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CHAPTER  1 


INTRODUCTION 

Optimal  periodic  control  theory  is  an  interesting  and 
relatively  unexplored  branch  of  optimization  theory.  There  are 
many  applications  identified  in  the  literature  and  contributions  to 
the  theory  are  increasing.  In  the  first  section  of  this  chapter 
preliminary  remarks  provide  the  general  background  and  motivation 
for  this  research.  After  a  statement  of  the  problem  and  its  assump¬ 
tions,  previous  theoretical  development  is  briefly  reviewed.  Then 
the  objectives  of  this  research  and  its  scope  are  outlined.  A 
description  of  the  notation  conventions  employed  concludes  the 
chapter. 

1 . 1  Preliminary  Remarks 

An  optimal  control  problem,  in  general,  can  be  reduced  to 
a  non-linear  two-point  boundary  value  problem  in  terms  of  the  state 
variables  and  Lagrange  multipliers.  Except  for  relatively  few  cases, 
an  analytical  solution  to  the  boundary  value  problem  is  untractable. 
One  of  the  exceptions  is  the  class  of  problems  satisfied  by  the 
static  equilibrium  solution.  For  many  practical  applications, 
reasonable  models  of  plant  processes  are  approximately  constants. 
Consequently,  many  control  systems  are  operated  "optimally"  by 
regulating  them  as  near  as  possible  to  an  equilibrium  state. 


Improved  performance  frequently  may  be  obtained  by  allowing 
the  state  of  one  of  these  regulated  processes  to  vary  with  time. 
Notably  in  the  field  of  Chemical  Engineering,  there  are  many  appli¬ 
cations  studied  in  the  literature  [1,  2,  3]*  in  which  periodic  or 
cycling  operation  provides  improved  performance  over  steady  state 
operation.  This  was  the  motivation  for  much  of  the  early  theoretical 
work  in  optimal  periodic  processes  by  Horn  and  Lin  [4]. 

Another  particularly  interesting  example  is  the  result 
obtained  by  Speyer  £ 5 ]  in  a  controversy  [5,  6,  7,  8]  over  an  aerospace 
problem.  He  showed  that  the  static  cruise  point  is  non-minimizing 
with  respect  to  range  factor  (fuel  rate  per  range  rate)  for  several 
dynamic  models  of  aircraft.  He  also  showed  numerically  that  operating 
along  a  cyclic  but  non-minimizing  path  provided  better  performance 
than  operating  at  the  static  cruise  point.  He  was  not  able  to  show 
what  the  optimum  path  was,  however,  and  to  date,  this  problem  remains 
unsolved. 

Later  numerical  attempts  [9]  to  obtain  a  periodic  solution 
to  this  problem  urging  standard  optimization  techniques  such  as 
steepest  ascent  and  conjugate  gradient  methods  failed  to  converge. 
Comparable  results  obtained  in  the  initial  research  for  this  work 
underscored  the  complexity  of  the  particular  application  and  motivat¬ 
ed  refocussing  the  present  research  effort  on  a  more  fundamental 

♦Numbers  appearing  in  brackets  refer  to  references  listed  in  the 
bibliography. 
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and  elucidating  example.  A  specific  problem  has  been  constructed 
which  is  composed  of  very  simple  dynamics  and  a  performance  index 
allowing  optimal  periodic  motion. 

1.2  Statement  of  the  Problem 

The  optimal  periodic  control  problem  is  a  subclass  of 
optimum  control  problems  that  has  two  distinguishing  features. 

First,  the  measure  of  performance  (performance  index  or  cost)  is  an 
average  over  one  period  of  the  performance  criterion.  Second, 
the  periodicity  constraint  requires  the  initial  and  final  states 
to  be  coincident.  The  following  formulation  and  assumptions  define 
the  general  optimal  periodic  control  problem  studied  in  this  work. 

Find  the  period,  T,  and  the  control,  u,  which  minimize  the 
performance  index, 

J  =  |(  L(x,u)dt,  (1.1) 

subject  to  the  constraints  of  the  system,  defined  by 

x  =  f (x , u)  ,  (1.2) 

and  to  the  periodicity  condition  applied  at  the  boundary, 

x(T)  =  x(0).  (1.3) 

The  assumptions  applicable  to  this  problem  are: 

(1)  The  period,  Te(0,+oo); 
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(2)  The  state  vector,  x(t)€XCRn,  where  X  is  the  set  of 
all  permissible,  continuous,  T-periodic  state  functions  in  the  n  dim¬ 
ensional  state  space; 

(3)  The  control  vector,  u(t)€UCRm,  where  U  is  the  set  of 
all  permissible,  piecewise  continous,  T-periodic  control  functions 

in  the  m  dimensional  input  function  space;  and 

(4)  The  functionals  f(-,*)  and  L(*,*)  and  their  first  two 
derivatives  with  respect  to  both  arguments  are  continuous.  Note  that 
the  system,  f,  and  the  performance  criterion,  L,  are  both  time  in¬ 
variant. 

1 . 3  Historical  Review 

One  of  the  earliest  contributions  to  the  theoretical 
development  of  optimal  periodic  control  is  the  work  of  Horn  and  Lin 
[4].  They  examine  chemical  engineering  steady-state  and  batch 
processes  that  can  be  improved  by  cycling.  Using  a  formulation  for 
the  problem  similar  to  the  previous  section,  they  derive  the  first 
order  necessary  conditions  for  optimality.  For  the  optimal  period, 
they  obtain  the  special  form  of  the  transversality  condition  which 
states  that  the  performance  index  must  equal  the  variational 
Hamiltonian  evaluated  along  the  extremal  path.  This  special  condition 
for  optimal  periodic  processes  has  not  been  previously  exploited  in 
numerical  procedures.  Another  early  contributor  to  this  theory  is 
Fjeld  [10]. 

The  practical  problem  of  establishing  whether  the  optimal 
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steady  state  process  can  be  improved  by  a  suitable  periodic  process 
has  been  analyzed  by  Bittanti,  Fronza,  and  Guardabassi  [11].  They 
establish  a  new  local  sufficiency  condition  which  involves  trans¬ 
forming  a  second  variation  quadratic  expression  into  the  frequency 
domain  and  examining  the  result  for  non-negativity  over  the  range 
of  all  positive  frequencies.  This  frequency  test  along  with  other 
equivalent  tests  are  presented  in  the  work  of  Willems  [12].  He 
intensively  treats  the  optimal  control  of  linear  systems  with  respect 
to  quadratic  but  not  necessarily  convex  performance  criteria. 

The  dynamic  programming  sufficiency  conditions  for  the 
optimality  of  periodic  control  processes  are  developed  in  a  later 
work  by  Maffezzoni  [13].  Global  conditions  are  derived  for  the 
optimal  control  from  the  Hamilton-Jacobi  equation.  Local  suffici¬ 
ency  condi tons  for  a  weak  minimum  are  derived  by  Bittanti,  Locatelli, 
and  Maffezzoni  [14]  for  periodic  processes  with  a  given  period.  They 
use  calculus  of  variations  techniques  similar  to  that  of  Bryson  and  Ho 
[15].  They  also  develop  sufficient  conditions  for  the  existence  of 
periodic  solutions  to  the  Riccati-type  equations  that  result  from 
the  second  variation  analysis. 

In  more  recent  work,  first  order  necessary  conditions  are 
reviewed  by  Gilbert  [16,  17]  for  a  very  general  formulation  of  the 
optimal  periodic  control  problem.  He  develops  relationships  between 
elements  of  the  solution  sets  for  the  optimal  periodic  control  problem 
and  optimal  steady-state  problem.  Bernstein  and  Gilbert  [18,  19] 
identify  an  overlooked  normality  condition  in  the  development  of 


the  frequency  test  [11].  They  also  expand  the  test  to  include  the 
more  general  formulation  of  the  problem  above. 
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Much  of  the  literature  on  optimal  periodic  control  is 
reviewed  by  Guardabassi,  Locatelli,  and  Rinaldi  [20]  in  an  important 
survey  paper.  Significant  open  aspects  of  the  problem  are  identified, 
some  of  which  are  specifically  addressed  in  this  work.  Another  review 
by  Bailey  [21]  emphasizes  applications  of  the  theory  to  a  particular 
chemical  engineering  problem. 

Important  theoretical  work  from  other  disciplines  should 
also  be  mentioned.  There  is  considerable  research  devoted  to  the 
study  of  periodic  solutions  (orbits)  for  dynamic  systems.  They 
date  back  to  at  least  1890  and  the  extensive  work  by  Poincare  [22] 
in  celestial  mechanics.  With  regard  to  finding  periodic  orbits, 
the  most  studied  dynamic  system  is  the  restricted  problem  of  three 
bodies.  An  excellent  and  thorough  review  of  important  work  on  this 
problem  is  provided  by  Szebehely  [23]. 

Since  his  earlier  work  [24]  on  periodic  orbits  near 
equilibrium  points  of  the  restricted  problem  of  three  bodies,  Broucke 
has  meticulously  studied  periodic  orbits  of  a  large  variety  of  fourth 
order  dynamic  systems.  Much  of  this  work  is  in  preparation  for 
publication.  Two  extensive  studies  of  periodic  solutions  to  other 
fourth  order  systems  are  by  Contopoulos  [25,  26]  and  Henon  [27]. 

A  common  element  in  all  of  these  studies  is  the  complex¬ 
ity  of  the  resulting  periodic  solutions.  Families  of  solutions 
exist  which  intersect  at  common  critical  solutions  (bifurcation 


points),  and  the  solutions  densely  pack  large  regions  of  state 
space.  Quasi-periodic  solutions  are  also  studied.  This  whole  area 
of  dynamics  is  currently  engrossed  in  intense  research  efforts. 

The  text  by  Meirovitch  [28]  provides  numerous  special 
properties  of  Hamiltonian  systems  with  periodic  coefficients. 

Properties  of  the  monodromy  matrix,  the  system  transition  matrix 
evaluated  over  one  period,  are  of  particular  interest  to  the 
present  work. 

1.4  Purpose  and  Scope 

The  main  objective  of  this  research  is  to  develop  a 
useful  understanding  of  the  structure  and  characteristics  of 
periodic  processes  resulting  from  the  optimal  periodic  control 
problem.  An  attempt  has  been  made  in  this  work  to  integrate  applicable 
portions  of  the  theory  and  related  experience  developed  in  the  fields 
of  analytical  dynamics  and  celestial  mechanics  with  that  existing  for 
this  problem  derived  under  the  discipline  of  optimal  control  theory. 

It  is  hoped  that  this  expansion  of  the  theory  will  lead  to  its  increas¬ 
ed  application  and,  in  particular,  to  the  development  of  improved  num¬ 
erical  techniques  for  analyzing  more  complex  systems. 

There  are  three  basic  parts  that  comprise  this  study. 


The  first  part  involves  the  general  theoretical  background  of  the 
problem.  Necessary  conditions,  sufficient  conditions,  and  tests 
for  the  optimality  of  continuous,  time-invariant,  periodic,  control 


systems  are  presented.  Special  properties,  studied  by  analytical 
dynamicists,  of  time-invariant,  Hamiltonian  systems  with  periodic 
constraints  are  reviewed.  Of  particular  importance  are  the  proper¬ 
ties  of  the  monodromy  matrix,  a  constant  matrix  determined  by  the 
transition  matrix  for  a  periodic  Hamiltonian  system  evaluated  over 
one  period.  Concluding  the  first  part,  contributions  to  the  general 
theory  developed  during  the  present  research  are  derived.  Included 
is  a  new  algebraic  Riccati  equation.  It  determines  the  initial 
values  for  periodic  solutions  to  the  Riccati  differential  equation 
generated  from  a  second  variation  analysis.  A  similarity  transform¬ 
ation  of  the  monodromy  matrix  separates  the  system  eigenvalues  into 
the  two  diagonal  blocks  of  the  resulting  matrix.  This  results  in 
the  necessary  condition  for  periodic  extremal  solutions  that  cor¬ 
responding  monodromy  matrices  have  no  distinct  eigenvalues  on  the 
unit  circle.  Local  sufficiency  conditions  are  then  extended  to  the 
free  period  case,  and  finally,  properties  of  a  periodic  regulator 
are  developed. 

In  the  second  part  of  the  study,  a  particular  control  problem 
is  defined  that  assures  the  existence  of  optimal  periodic  solutions  for 
some  range  of  values  of  a  parameter  weighting  the  control.  A 
fairly  comprehensive  numerical  investigation  of  this  control  problem 
is  conducted.  Periodic  solutions  to  the  Euler-Lagrange  equations 
resulting  from  the  first  order  necessary  conditions  are  computed. 

They  form  one  parameter  families  which  intersect  at  common  "critical" 
solutions  or  bifurcation  points.  Several  solutions  are  examined 


along  and  between  families  to  identify  characteristic  differences 
that  distinguish  the  families.  The  necessary  condition  for  an  optimal 
period  is  applied  which  is  satisfied  by  many  of  the  solutions. 
Sufficiency  conditions  are  checked  for  these  solutions  by  integrat¬ 
ing  the  corresponding  matrix  Riccati  equation  over  one  period  to 
determine  whether  the  Riccati  variable  exists  during  the  entire 
interval.  Finally,  a  neighboring  optimum  feedback  control  system 
is  developed  to  operate  near  the  optimal  periodic  process.  The 
initial  state  is  perturbed,  demonstrating  the  behavior  of  the 
periodic  regulator. 

The  final  part  of  the  study  is  the  development  of  an 
approximate  analytical  solution  for  the  particular  control  problem 
previously  defined.  A  perturbation  technique  is  used  that  results 
in  a  solution  in  the  form  of  an  asymptotic  series  expansion.  The 
necessary  condition  for  optimal  period  is  applied  leading  to  expres¬ 
sions  in  the  perturbation  parameter  and  the  time  for  the  optimal 
path,  period,  control,  and  performance  index.  These  results  are 
compared  to  a  linear  analysis  near  the  static  solution  and  to 
earlier  numerical  results  obtained  at  the  minimum  solution. 

Conclusions  and  recommendations  for  further  research 
complete  the  study. 

1 .5  Notation  Conventions 

Set  theory  notation  is  used  sparingly  in  this  study. 

Symbols  used  and  their  intended  definitions  are  listed  below: 
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(1)  e  signifies  "is  an  element  of", 

(2)  c  signifies  "is  a  subset  of",  and 

(3)  Rn  signifies  a  real,  Euclidean  n-dimensional  hyperspace. 
Thus  the  notation  T€(0,+«)  means  that  T  is  an  element  of  the  open  set 
0  <  T  <  +»,  and  the  notation  x(*)SXCRn  means  that  x  is  an  n-vector 
of  the  set  of  n-vectors  X  which  is  a  subset  of  all  vectors  defined  in 
n-dimensional,  real,  Euclidean  hyperspace. 

Vector  and  matrix  notation  are  used  throughout  this  study. 

The  following  conventions  are  used: 

(1)  Vectors  are  generally  indicated  by  lower  case  letters 
and,  unless  otherwise  indicated,  are  column  vectors.  Components  of 
vectors  are  denoted  by  subscripts.  Thus,  the  n-vector  x  and  its  com¬ 
ponents  are: 


(2)  Matrices  are  indicated  by  upper  case  letters  and 
sometimes  by  elements  enclosed  in  brackets .  Elements  of  matrices 
are  denoted  by  double  subscripts.  Thus,  the  n  row  by  m  column 
matrix  A  and  its  elements  are  denoted  by 

«  •••  3 . 

11  lm 

a  ....  a 
nl  am 


A 


(3)  Scalars  are  easily  recognized  as  such,  otherwise 


they  are  specifically  identified  when  used.  Examples  are  the 
independent  variables,  t  and  X,  for  time  and  extended  time,  and 
the  indexing  variables  i,  j,  and  k.  Others  are  the  performance 
index,  J,  and  the  variational  Hamiltonian,  H. 

(4)  The  transpose  of  a  vector  or  matrix  is  denoted  by  the 

T 

superscript  T,  thus  A  denotes  the  transpose  of  matrix  A. 

(5)  The  inverse  of  a  square  matrix  is  denoted  by  a  super¬ 
script  -1,  thus  B  *  denotes  the  inverse  of  matrix  B. 

(6)  The  first  partial  derivative  of  a  scalar  with  respect  to 
a  vector  is  a  row  vector  designated  by 

H=3H=r§2<M  3H] 

x  ~  3x  '  [  3Xj  8X2’  ’  *  3xn  J  ' 


(7)  The  second  partial  derivative  of  a  scalar  with  respect 
to  a  vector  is  a  matrix  designated  by 


H 


xu 


Note  that  the  transpose  of  the  matrix  H  is  H 

r  xu  ux 

(8)  The  first  partial  derivative  of  a  vector  with  respect  to 
another  vector  is  a  matrix  designated  by 
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f 

x 


3fl  .  .  . 

8fi 

8xl 

dx 

a'f 

n  .  .  . 

8f 

n 

axi 

3x 

n 

(9)  The  first  partial  derivative  of  a  matrix  with  respect 
to  a  vector  is  designated  by 


f 

xu 


2 

a^f . 

_ i 

i 

3u.  8x . 
k  J 


where  i  and  j 


1, 


n  and  k  =  1 ,  . . . ,  m. 


below: 


Other  notation  conventions  used  in  this  study  are  listed 


(1)  The  dot  convention  is  used  to  designate  a  total  deriv¬ 
ative  with  respect  to  the  independent  variable,  generally  time.  Thus 
the  following  notation  is  equivalent, 


x 


dx 
dt  ‘ 


(2)  The  derivative  operator,  D  ,  is  used  in  chapter  six  to 
designate  the  n^  derivative  with  respect  to  the  expanded  time,  T, 


2 

d  x  _  _2 

— j  =  D  x. 

dt 


(3)  Various  notation  conventions  are  used  in  the  literature 
to  identify  the  variations  of  a  function  or  functional.  In  this  study 
only  the  variation  operators  d  and  6  are  used,  conforming  to  the 
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predominant  usage  present  in  the  literature.  Unfortunately  d  is 
also  used  to  denote  a  differential  along  the  path  which  provides  a 
source  of  confusion  in  the  literature.  The  distinction  between  the 
variation  operators  is  that  d  represents  a  total  variation  while  6 
represents  a  variation  with  time  held  constant.  The  following  expres¬ 
sions  provide  the  functional  relationship  between  them;  figure  1.1 
graphically  identifies  this  relationship, 

dx(tf)  =  x(tf+dtf)  -  x°(tf), 

=  x(tf)  +  x(tf)dtf  +  |x(tf)dt^+. . .-  x°(tf ) , 

=  6x(t^)  +  x°(t^)dt^  +  higher  orders.  (1.4) 

In  equation  (1.4)  and  figure  1.1,  x°  is  the  extreme  (reference)  tra¬ 
jectory,  x  is  a  neighboring  (comparison)  path,  and  x°  is  the  result, 
to  first  order,  of  expanding  x  in  a  Taylor  series  about  x°.  In  later 
chapters  the  superscript  "o"  designates  the  optimal  trajectory  or 
optimal  solution. 


*  " 


CHAPTER  2 


CONDITIONS  FOR  OPTIMALITY 


The  conditions  for  optimality  of  a  control  process  may  be 
derived  by  at  least  two  different  methods.  The  classical  calculus  of 
variations  approach  is  used  in  many  texts  [15,  29,  30,  31]  to  derive 
the  conditions  for  local  optima;  whereas,  Bellman's  dynamic  program¬ 
ming  technique  [32,  33,  34],  using  the  Hamilton-Jacobi  theory,  also 
provides  global  results.  In  the  first  section  of  this  chapter  the 
existing  first  order  necessary  conditions  for  optimal  periodic  pro¬ 
cesses  are  derived  using  the  calculus  of  variations.  The  next  two 
sections  review  existing  sufficient  conditions  for  optimality  derived 
from  each  method.  The  frequency  test,  which  determines  whether  or 
not  a  periodic  process  can  improve  steady-state  performance  in  a 
control  system,  is  presented  in  the  final  section. 


2.1  First  Order  Necessary  Conditons 


The  optimal  periodic  control  problem  defined  in  the  first 
chapter  by  equations  (1.1)  through  (1.3)  is  rewritten  here  for  con¬ 
venience.  Find  the  period,  T,  and  the  control,  u,  which  minimize 
the  performance  index, 


J  "  T 


fT 

1  L(x,u)dt, 

Jo 


subject  to  the  following  system  constraint, 
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*  =  f(x,u),  (1.2) 

and  the  periodicity  condition  applied  at  the  boundary, 

x(T)  =  x(0) .  (1.3) 

The  system  constraint  and  periodic  boundary  conditions 
can  be  adjoined  to  the  performance  index  with  Lagrange  multipliers 
forming  the  augmented  performance  index, 

5  s  vT<P(x(T),x(0))  +  \  j  [H(x,u,A)-XTx]dt.  (2.1) 

The  new  functions  introduced  by  equation  (2.1)  are  the  variational 
Hamiltonian, 

H(x,u,A)  s  L(x,u)  +  \Tf(x,u),  (2.2) 

and  the  boundary  condition  null  identity, 

4*(x(T)  ,x(0) )  =  x(T)  -  x(0).  (2.3) 

The  n-vector  Lagrange  multipliers,  V  and  X,  multiply  the  null 
identities  formed  from  the  boundary  condition  equation  (1.3)  and 
the  system  constraint  equation  (1.2)  respectively. 

First  order  necessary  conditions  for  the  optimal  periodic 
control  are  now  derived.  From  fundamental  arguments  in  the  calculus  of 
variations  the  first  variation  of  the  augmented  performance  index  must 
vanish.  Written  as  the  sum  of  three  terms,  the  first  variation  is 
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dJ  =  d(vTY)  +  d(i)  1  [H-XTx]dt  + 


[H-XTx]dt  +  |d(  J  (H-XTx 


x]dt) .  (2.4) 


Expanding  the  initial  term  of  the  first  variation  provides 

d(vT40  =  dvT¥  +  vVx(T)dx(T)  +  Vx(0)dx(0)]. 

This  may  be  simplified  by  recognizing,  from  equation  (2.3),  that 
1*x^j  =  1  and  =  “1-  Substituting  these  into  the  above  expres¬ 

sion  for  the  first  term  gives 


d(vT40  =  dvTY  +  vT[dx(T)  -  dx(0) ) . 


(2.5) 


Operating  as  indicated  on  the  second  term  of  equation  (2.4)  yields 


dd)(  [H-XTx]dt  =  [H-XTx]dt  . 

1  JO  T  JO 


Simplifying  this  by  using  equations  (1.1),  (1.2),  and  (2.2)  gives 


d(|)J  [H-XTx]dt  =  -  JpJ. 


(2.6) 


The  expansion  of  the  last  term  of  equation  (2.4)  is  common  to  free 
terminal  time  optimal  control  problems, 


|d(  j  V-XTx 


x]dt)  =  |(H-XTx)dT 


•ti 


(Hx6x  +  Hu6u  +  6X^(H^-x)  -  X^6x]dt. 
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Integrating  the  last  tern  of  the  right  hand  integral  by  parts  and  using 
the  relationship  (valid  to  first  order),  6x  =  dx  -  xdt,  reduces  the 
previous  equation  to 


[H-XTx]dt) 

0 


^HdT  +  |(XTdx) 


t  =  T 
t  =  0 


+ 


[(Hx+XT)6x  +  Hu6u  +  6XT(H^-x) ]dt  . 


(2.7) 


Combining  the  results  for  the  three  terms  given  by  equations  (2.5), 
(2.6),  and  (2.7)  provides  the  following  expression  for  the  first 
variation, 


dJ 


dvTY  +  |(H-J)dT  +  [(vT+  |\T(t))dx(t)J 


t  =  T 

t  =  0 


+ 


[(H  +XT)6x  +  H  6u  +  6XT(HT-x)]dt. 
0  x  u  A 


(2.8) 


The  criteria  that  the  first  variation  must  be  zero  when 
evaluated  along  an  extremal  path  will  now  be  examined.  The  terms 
which  include  the  Lagrange  multiplier  variations  vanish  since  their 
coefficients  are  the  null  identities  previously  defined.  Recognizing 
that  using  Lagrange  multipliers  permits  treating  all  other  variables 
as  though  they  were  independent  allows  the  remaining  coefficients  of 
variations  in  equation  (2.8)  to  be  isolated  by  appropriately  selecting 


admissable  comparison  paths.  Thus  each  of  these  coefficients  must  also 
equate  to  zero  when  evaluated  along  an  extremal  path. 

The  Euler-Lagrange  equations  result  from  requiring  the 
coefficients  of  the  variations  in  the  integral  of  equation  (2.8) 
to  vanish.  They  are: 


x  =  H* 

=  f(x,u). 

(2.9) 

A  =  -HT 

X 

=  -f*(x,u)X  -  l7(x,u), 

(2.10) 

T 

0  =  H  = 
u 

f^(x,u)A  +  L*(x,u). 

(2.11) 

The  fixed  point  relationships  are  established  by  equating  the  remain¬ 
ing  coefficients  of  the  equation  to  zero.  The  coefficient  in  the 
first  term  is  the  boundary  null  identity  and  gives  back  the  pre¬ 
scribed  boundary  conditions, 

x(T)  =  x(0).  (2.12) 

The  natural  boundary  conditions,  i.e.,  transversality  condition,  comes 
from  the  remaining  terms, 

\(T)  =  A(0)  =  -Tv,  (2.13) 

H(T)  =  H(0)  =  J°.  (2.14) 

Any  periodic  solution  to  the  two  point  boundary  value 
problem,  equations  (2.9)  through  (2.13),  is  an  extremum  of  the  problem 


when  the  period  is  specified.  The  condition  (2.14)  relating  the 
variational  Hamiltonian  and  the  performance  index,  evaluated  along 
the  optimal  path  is  the  special  condition,  first  derived  by  Horn  and 
Lin  [4] ,  for  testing  the  optimal  period.  For  the  case  where  H  is  auto¬ 
nomous,  which  is  assumed  here,  the  Hamiltonian  is  a  constant  on  any 
extremal  path,  and  equation  (2.14)  is  valid  for  all  times.  This 
condition  is  obviously  not  affected  by  the  particular  choice  of  the 
initial  time. 

2.2  Sufficient  Conditions  for  a  Local  Minimum 

A  set  of  sufficient  conditions  for  local  optimality  of  the 
periodic  control  problem,  defined  by  equations  (1.1),  (1.2),  and  (1.3) 
with  the  period  T  given,  is  reviewed.  These  conditions  were  initially 
derived  by  Bittanti,  Locatelli,  and  Maffezzoni  [14]  using  first  and 
second  variations  and  properties  of  periodic  functions.  In  the  fol¬ 
lowing  restatement  of  the  sufficient  conditions,  the  function  4>A(T,0) 
is  the  transition  matrix  for  the  system  y  =  Ay  evaluated  over  one 
period,  and  the  control  u°(0  is  a  piecewise-continuous ,  T-periodic 
function  in  the  input  function  space.  The  superscript  "o"  designates 
the  optimum  solution. 

The  control  u°(»)  is  a  local  minimum  if  the  following  con¬ 
ditions  are  satisfied: 

(1)  The  Euler-Lagrange  equations,  (2.9)  through  (2.11), 
the  prescribed  boundary  condition,  (2.12),  and  the 
transversality  condition,  (2.13),  are  satisfied; 
u°(*)  is  a  solution  to  (2.11)  for  0  2  t  2  T,  and 


x°(*)  and  X°(*)  are  the  corresponding  solutions  to 
the  resulting  two-point  boundary  value  problem; 
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(2)  No  eigenvalue  of  <!>  o(T,0)  is  equal  to  one; 

^x 

(3)  The  strong  form  of  the  Legendre  condition  is  satisfied, 


H  >  0,  for  0  <  t  <  T; 
uu  -  - 


(4)  A  bounded  symmetric  solution  to  the  following  Riccati- 


type  equation  exists  for  0  <  t  <  T, 


P  =  -  PA  -  A  P  +  PBP  +  C, 


(2.15) 


where  A  = 


f°  -  f°(H°  )-1H°  , 
x  u  uu  ux’ 


f°(H°  )”1f°T  and 
u  uu  u’ 


-H°  +  H°T(H°  )-1H°  , 

XX  ux  uu  ux’ 


subject  to  the  periodic  boundary  condition 


P(T)  =  P(0) ;  and 


(5)  No  eigenvalue  of  ^(TiO)  is  equal  to  one,  where 


Z  =  f°  -  f°(H°  )_1(H0  +  f°TP). 

X  u  uu  ux  u 


(2.16) 


(2.17) 


It  will  be  shown  in  chapter  four  that  the  last  condition  (5) 
as  stated  in  [14],  is  never  true  since  <t>^(T,0)  always  has  a  unit  eigen¬ 
value.  A  necessary  condition  for  optimality  is  derived  from  the 


t 
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eigenvalues  of  this  matrix.  Also,  an  algebraic  Riccati  equation  is 
determined  that  provides  the  initial  conditions  for  periodic  solutions 
to  equation  (2.15)  in  condition  (4). 

2.3  Hamilton-Jacobi  Theory 

By  extending  the  Hamilton-Jacobi  theory  to  periodic  optimiza 
tion  problems,  Maffezzoni  [13]  derived  a  sufficient  condition  for  opti 
mality  of  the  problem  defined  by  equations  (1.1),  (1.2),  and  (1.3) 
with  the  period,  T,  free.  It  is  assumed  that  the  solution  V(t,x)  to 
the  Hamilton-Jacobi  equation, 

Vt(t,x)  +  H(x,Vx(t,x),u)  =  0,  (2.18) 

is  twice  continuously  diffentiable  with  respect  to  both  of  its  argu¬ 
ments  throughout  phase  space.  It  is  also  assumed  that  the  varia¬ 
tional  Hamiltonian  is  regular,  i.e.,  that  it  has  a  unique  minimum 
with  respect  to  the  control  for  given  values  of  its  arguments.  With 
these  assumptions,  a  statement  of  the  sufficient  condition  follows. 

If  a  real  function  C(*)  exists  relating  solutions  to  (2.18) 
in  the  form 


V(0,x)  -  V(T,x)  =  C(T) ,  (2.19) 

ic 

and  u  (x,Vx(t,x))  is  the  control  that  minimizes  the  Hamiltonian  in 
equation  (2  18),  for  all  x€X,  then  the  optimal  control  relative  to  x 
and  the  specified  T  is 


u°(t)  =  u*(x°(t),Vx(x°(t),t)), 


(2.20) 


<  V.  J  v  -  *  -  - 
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where  u°  and  x°  denote  the  optimal  control  and  state  relative  to  the 
period . 

A  corollary  relationship  to  the  previous  condition  is  also 
developed  in  this  work.  Under  the  previous  assumptions  the  optimal 
performance  relative  to  x  and  T  can  be  expressed  as 

J°  =  V<0,«V  =  Cp.  (2.21) 

2.4  The  Frequency  Test 

The  question  of  whether  the  optimal  steady-state  operation 
of  a  given  plant  can  be  improved  by  a  cycling  or  periodic  process  is 
an  important  preliminary  concern.  This  question  has  been  considered 
in  much  of  the  previous  work  on  optimal  periodic  control  processes. 

From  a  relaxed  steady-state  control  theory  approach,  viola¬ 
tion  of  the  maximum  principle  by  the  optimal  steady  state  control 
leads  to  an  optimal  solution  of  the  bang-bang  or  chattering  type. 
Reference  should  be  made  to  Horn  and  Lin  [4],  Bailey  and  Horn  [1], 
and  Gilbert  (16,  35]  for  further  discussion  and  examples. 

A  second  variation  analysis  by  Bittanti,  Fronza,  and 
Guardabassi  [11]  established  a  frequency  domain  local  sufficiency 
condition  for  determining  whether  the  steady  state  control  of  a  system 
can  be  improved  by  cycling.  However,  the  condition  does  not  provide 
sufficient  information  to  determine  what  cyclic  or  periodic  process 
is  best.  An  overlooked  normality  (controllability)  condition  was 
identified  by  Bernstein  and  Gilbert  [18]  in  a  somewhat  more  general 
treatment  of  the  same  problem. 
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The  frequency  test  is  stated  here  without  proof.  For  the 
optimal  periodic  control  problem  defined  by  equations  (1.1),  (1.2), 
and  (1.3),  where  the  system  is  assumed  to  be  controllable,  let 

A  h  fx(x°,u°),  B  =  fu(x°,u°), 


P  s  HXX(X°’U°  x0)’ 


Q  =  Hxu(x°,u°,A°), 


R 


H  (x°,u°,A°) 
uu  >  y  , 


(2.22) 


and  no  eigenvalue  of  A  has  a  zero  real  part.  Willems  allows  any  matrix 
A  in  his  work  [12].  The  superscript  on  x,  u,  and  X  designates  the 
optimal  solution  in  the  class  of  constant  solutions  of  the  optimal 
periodic  control  problem.  Form  the  (n  x  n)-Hermitian  matrix, 

f7(u0  =  GT(-jtu)PG(juj)  +  QTG(jw)  +  GT(-jw)Q  +  R,  (2.23) 


where  G(s)  =  (sI-A)  ^B,  ut  is  frequency,  and  j  =  V"l-  If  n(u>)  is  not 
negative  definite  for  all  ui,  then  the  optimal  control  for  the  problem 
is  not  constant,  but  belongs  to  the  class  of  time  varying  periodic 
controls . 
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CHAPTER  3 


HAMILTONIAN  SYSTEMS  WITH  PERIODIC  SOLUTIONS 

The  optimal  periodic  control  problem  can  be  expressed  in 
the  standard  canonical  form  of  a  Hamiltonian  system  by  functionally 
eliminating  the  control  from  the  Euler-Langrange  equations.  The 
transition  matrix  of  this  Hamiltonian  system  evaluated  over  one  or 
an  integral  multiple  of  the  period  and  the  corresponding  eigenvalues 
of  this  matrix  have  properties  which  are  exploited  in  the  theoretical 
developments  of  the  next  chapters.  These  important  properties,  as 
well  as  the  fundamental  properties  of  the  Hamiltonian  systems  from 
which  they  are  derived,  are  reviewed  in  this  chapter.  There  are  a 
vast  number  of  references  available  that  address  different  portions 
of  the  material  covered  here  in  various  degrees  of  completeness. 

The  intent  of  this  chapter  is  to  include  in  one  place,  in  a  consistent 
notation,  a  statement  and  derivation  of  properties  that  are  fundamental 
to  new  theory  developed  in  later  chapters  of  this  work. 

3. 1  Hamiltonian  Structure  of  the  Problem 

In  general  the  optimal  periodic  control  problem  can  be 
expressed  in  the  familiar  canonical  form  of  Hamilton's  equations  of 
motion.  The  necessary  condition  for  all  extrema,  that  the  first 
variation  of  the  performance  index  vanish,  yields  the  Euler-Lagrange 
equations , 
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H^(x,u,A), 

(3.1) 

-H*(x,u,A), 

(3.2) 

H^(x,u,A), 

(3.3) 

subject  to  the  periodicity  conditions  at  the  boundary,  where  T  is  the 
period, 


x(T)  =  x(0)  and  A(T)  =  A(0) .  (3.4) 

An  expression  for  the  control,  u,  in  terms  of  the  state 
variables  and  Lagrange  multipliers,  x  and  A,  may  be  determined  (at 
least  implicitly)  from  the  last  of  the  Euler-Lagrange  equations  (3.3). 
Eliminating  the  control  from  the  remaining  two  equations  reduces  the 
system  of  equations  to  the  desired  canonical  form, 


H^(x,u(x,A),A), 

(3.5) 

-H^(x,u(x,A),A). 

(3.6) 

In  terms  of  the  dynamicist,  the  variational  Hamiltonian 
H(x,A)  and  the  two  n-vectors  of  the  state,  x,  and  the  Lagrange  multi¬ 
pliers,  A,  are  respectively  the  Hamiltonian  function  and  the  two  phase 
space  n-vectors  of  coordinates  and  momenta. 

It  is  convenient  to  make  a  further  simplification.  First 


define  the  (2n  x  2n)  skew  symmetric  matrix,  K,  which  is  commonly 
called  the  fundamental  symplectic  matrix, 
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(3.7) 


where  I  is  the  n  square  identity  matrix.  By  inspection,  the  funda¬ 
mental  symplectic  matrix  satisfies 

K~ 1  =  KT  =  -K.  (3.8) 

Now  Hamilton's  equations  may  be  written  in  the  simplified  form 

Y  =  KHj  =  Y(y),  (3.9) 

where  y  is  the  2n  vector  composed  of  the  phase  space  components  x 
and  A  and  is  subject  to  the  periodicity  condition 

y(T)  =  y(0) .  (3.10) 

From  equation  (3.9)  it  is  apparent  that  the  tangent  to  a  phase  space 
trajectory  is  a  function  of  the  position  vector,  y,  only.  Hence,  there 
is  only  one  trajectory  through  a  point  in  phase  space,  and  each  trajec¬ 
tory  is  fixed.  It  will  be  convenient  later  to  identify  these  trajec¬ 
tories  (periodic  solutions  to  the  Euler-Lagrange  equations)  by  their 
initial  conditions.  Note  that  the  system  is  assumed  autonomous;  that 
is,  the  Hamiltonian  is  not  an  explicit  function  of  the  time. 

Another  well  known  property  of  an  autonomous  system  is  that 
the  Hamiltonian  is  a  constant  of  the  motion.  This  is  easily  shown  by 
substituting  in  the  expression  for  the  time  derivative  of  the 
Hamiltonian, 
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H(x,u,\)  =  H  x  +  Hi  +  H  u, 

x  A  u 

the  equations  of  motion,  (3.1)  through  (3.3),  which  gives 


(3.11) 


H(x,A)  =  HH*  +  H^(-H*)  .  (3.12) 

Since  the  scalar  inner  product  is  invariant  to  the  transpose,  equation 
(3.12)  reduces  to 


(3.13) 


and  the  Hamiltonian  is  a  constant  of  the  motion. 

The  first  integral  may  be  used  to  reduce  by  one  the  order 
of  the  differential  equations  that  define  the  system.  Any  variables 
in  the  Hamiltonian  function  is  expressable  in  terms  of  a  constant 
parameter  and  the  remaining  variables.  Given  any  set  of  initial 
conditions  the  constant  parameter  is  uniquely  determined.  A  further 
reduction  in  the  order  of  the  system  is  possible  by  eliminating  the 
independent  variable,  time,  from  the  equations.  The  computational 
scheme  employed  in  this  work  to  determine  periodic  solutions  to  the 
Euler-Lagrange  equations  exploits  the  first  of  these  relationships  and 
in  effect  reduces  the  degree  of  the  problem  by  one  order. 

3.2  Symplectic  Properties  and  the  Transition  Matrix 

Some  fundamental  characteristics  of  all  Hamiltonian  systems 
are  associated  with  the  properties  of  symplectic  matrices.  To  begin 
this  section,  the  definitions  of  two  types  of  synxplectic  matrices  and 


their  corresponding  properties  are  provided.  Then  the  transition 
matrix  of  a  Hamiltonian  system,  along  with  some  of  its  basic  proper¬ 
ties,  are  derived.  Finally,  it  is  shown  that  the  transition  matrix 
of  a  Hamiltonian  system  is  symplectic. 

First,  the  (2n  x  2n)  matrix,  B,  is  defined  to  be  a  symplec¬ 
tic  matrix  if 

BKBT  =  K,  (3.14) 

where  K  is  the  fundamental  symplectic  matrix  (3.7)  defined  in  the 
previous  section.  Second,  the  (2n  x  2n)  matrix,  C,  is  a  skew  symplec¬ 
tic  matrix  if 

KCK-1  =  -CT.  (3.15) 


With  the  use  of  the  symplectic  matrix  relationships  expressed  in  (3.8) 
many  equivalent  forms  of  (3.14)  and  (3.15)  may  be  generated. 

To  illustrate  properties  associated  with  these  matrices, 
partition  each  matrix,  B  and  C,  into  four  (n  x  n)  submatrices  which 
then  can  be  expressed  as 


‘Bll 

B12_ 

and  C  = 

'Cu 

Cl2‘ 

_B21 

®22 

C21 

C22_ 

(3.16) 


Substituting  the  partitioned  matrices  (3.16)  into  equations  (3.14)  and 
(3.15)  provides  the  following  properties  for  the  symplectic  matrix  B, 


*  *  •* 
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B22BIl  ’  <B22B21)T. 


B1IB22 


B  BT  =  Ib 
B12B21 


and  for  the  skew  symplectic  matrix  C, 


C12  =  C12’ 
C21  *  C21’ 


(3.18) 

(3.19) 


(3.20) 

(3.21) 

(3.22) 


One  important  result  of  equation  (3.22)  is  that  the  trace  of  a  skew 
symplectic  matrix  is  zero. 

Returning  to  the  Hamiltonian  system,  first  write  its  linear 
variational  equations.  This  may  be  done  by  considering  the  result  of 
a  small  perturbation,  6y,  from  a  periodic  solution,  y,  to  the  equations 
of  motion  (3.9).  Substituting  the  perturbed  solution  into  the  equa¬ 
tions  of  motion  gives 

y  +  6y  =  KHy(y+6y),  (3.23) 

which  may  be  expanded  in  a  Taylor  series  about  the  nominal  solution, 
y,  resulting  in 

y  +  8y  =  K[Hy  +  H^fiy  +  order(6y2)].  (3.24) 

Neglecting  6y  terms  of  order  two  and  higher  and  using  equation  (3.9) 
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in  the  result  gives  the  familiar  form  of  the  linear  variational 
equations , 


6y  =  A(t)6y, 


(3.25) 


where  A(t)  =  KH  and  the  matrix  of  second  partials,  H  ,  is  evaluated 

yy  F  ’  yy’ 

on  the  nominal  solution  y(t).  Note  that  the  matrix  H  is  a  real, 

yy 

symmetric  matrix  with  each  element  either  periodic  of  period  T  or 
constant.  Note  also  that  the  matrix  A(t)  is  skew  symplectic  as  may  be 
shown  by  substituting  the  matrix  into  equation  (3.15)  and  reducing  the 
result  to  the  identity  below: 


[  )K_1  = 

-(KH  ) 

yy 

yy 

T 

i 

« 

5* 

> 

II 

-H  K 

yy 

yy 

T 

-h  ir  = 

* 

i 

yy 

yy 

(3.26) 


The  relationships  for  K  expressed  in  equation  (3.8)  and  the  symmetry 
of  were  used  in  the  above  equations  (3.26). 

A  fundamental  solution  matrix,  <)>,  may  be  constructed  by  form¬ 
ing  a  square  matrix  composed  of  2n  columns  of  independent  solutions  to 
the  variational  equation  (3.25).  The  solutions  are  real  and  in  general 
not  periodic.  This  matrix  is  also  not  unique;  the  relationship  between 
any  two  fundamental  matrices  being 


♦  j(t)  =  02(t)C, 


(3.27) 


where  C  is  a  constant  matrix  for  all  t.  The  fundamental  matrix,  4>j, 


is  called  a  principal  fundamental  matrix  when  the  constant  matrix 


is  the  inverse  of  $2  evaluated  at  some  reference  time  tg.  This  very 
useful,  special  case  of  is  also  called  the  transition  matrix, 
typically  expressed  as  4»(t,tg).  Rewriting  equation  (3.27)  in  the 
following  form  defines  the  transition  matrix, 

4>(t,t0)  =  *(t)<f1(t0).  (3.28) 

Note  at  the  initial  time,  t  =  tg,  the  transition  matrix  reduces  to  the 
identity  matrix. 

From  its  construction  it  is  obvious  that  the  transition  matrix 
satisfies  the  variational  equation  (3.25),  thus 

<Kt,t0)  =  A(t)*(£,t0).  (3.29) 

Two  useful  rules  for  operations  on  transition  matrices  are  now  stated. 
The  composition  rule  shows  how  transition  matrices  may  be  combined, 

4>(t2,t0)  =  <Kt2,tj)  ♦(t1,tQ).  (3.30) 

The  inverse  of  a  transition  matrix,  is  functionally  obtained  by  inter¬ 
changing  its  arguments, 

*l>"1(t2,t1)  =  +(t1*t2)-  (3.31) 

Both  rules  are  easily  verified  by  application  of  the  definition  of  the 
transition  matrix  given  by  equation  (3.28), 

The  following  relationship  between  the  determinant  of  the 
transition  matrix  and  the  trace  of  the  coefficient  matrix  A(t)  is 
attributed  to  a  combination  of  Abel,  Jacobi,  and  Liouville  [36] 


33 


det  4>(t ,tQ) 


exp[ 


tr{A(o)}dcr] 

0 


(3.32) 


For  a  Hamiltonian  system  it  has  been  shown  in  (3.25)  and  (3.26)  that 
the  coefficient  matrix  A(t)  from  the  variational  equations  is  skew 
symplectic,  hence  the  trace  of  the  matrix  A(t)  is  zero.  From  equation 
(3.32),  the  determinant  of  the  transition  matrix  of  a  Hamiltonian 
system  is  one. 

An  important  property  of  the  transition  matrix  of  a 
Hamiltonian  system  is  that  it  is  symplectic.  As  a  consequence,  the 
transition  matrix,  (Kt.tg),  must  satisfy  equation  (3.14)  or  an  equiva¬ 
lent  form,  such  as 

KoV1  '  ♦"1.  (3.33) 

T  -1 

This  equality  may  be  verified  by  showing  that  K4»  K  is  a  solution  to 
the  adjoint  system, 

i'^t,^)  =  -<t>_1(t,t0)A,  (3.34) 


where  A  is  the  coefficient  matrix  from  (3.25).  The  adjoint  equations 
may  be  determined  by  operating  on  the  identity 

♦(t.t^V,^)  =  I.  (3.35) 

Taking  the  time  derivative  and  rearranging,  results  in  the  following 
relationship , 


(3.36) 


Using  equation  (3.29)  to  eliminate  ♦  in  the  left-hand  side  of  equation 
(3.36)  and  then  premultiplyiog  both  sides  of  the  equation  by  -♦  1  gives 
the  differential  equation  for  the  adjoint  system  previously  identified 
in  equation  (3.34). 

-1  T  -1 

To  verify  equation  (3.33),  replace  4>  with  K<fr  K  in  equa¬ 
tion  (3.34).  This  gives  the  following  result. 


•T  -1  T  -1 

K*XK  1  =  -K4»XK  A. 


(3.37) 


Eliminating  in  equation  (3.37)  using  equation  (3.29)  and  postmulti- 


plying  both  sides  of  the  resulting  equation  by  K  gives 


K<frTAT  =  -KOV1/^. 


(3.38) 


Using  the  skew  symplectic  property  of  A  given  by  equation  (3.15), 

T  -1 

the  above  equation  is  reduced  to  an  identity.  This  shows  K$  K 
satisfies  the  adjoint  equations. 

Investigating  initial  conditions,  when  t  =  tg,  the  follow¬ 
ing  relationships  of  the  transition  matrix  may  be  written, 


♦(W  =  *T<vv  -  ♦'W  =  (3-39) 


K^tg.t^K'1  =  KIK-1  =  I. 


(3.40) 


T  -i  -1 

Since  K$  K  and  ♦  satisfy  the  same  adjoint  equation  (3.34)  and 
they  are  equal  at  one  point  (3.40),  they  are  equal  everywhere.  Hence 


,  I  V/'-  V 


m 
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the  equality  of  equation  (3.33)  is  verified  and  the  symplectic  prop¬ 
erty  of  the  transition  matrix  of  a  Hamiltonian  system  has  been  demon¬ 
strated. 

3-3  Eigenvalues  of  the  Transition  Matrix 

The  basis  of  many  of  the  properties  of  a  linear  system  is 
its  eigenvalues.  They  provide  an  intrinsic  characterization  of  a 
matrix.  Of  principle  interest  are  the  relationships  associated  with 
a  special  matrix,  the  transition  matrix.  In  this  section  some  funda¬ 
mental  concepts  are  reviewed,  a  few  useful  relationships  between 
matrices  and  their  eigenvalues  are  presented,  and  finally  an  important 
property  of  the  eigenvalues  of  a  transition  matrix  is  established. 

Generally  the  eigenvalues  of  a  matrix  are  introduced  by  ref¬ 
erence  to  the  following  matrix  relationship, 

*£  =  MC,  (3.41) 

where  <t>  is  a  known  (2n  x  2n)  matrix,  £  is  a  nonzero  2n  vector,  and  p 
is  a  constant  scalar  parameter.  The  particular  values  of  p  that  satisfy 
this  relationship  are  the  eigenvalues  of  4>.  The  vectors  £  correspond¬ 
ing  to  each  of  the  eigenvalues  are  the  eigenvectors.  They  are  deter¬ 
mined  in  direction  only  and  form  a  set  of  basis  vectors  which  define 
the  2n  phase  space  of  the  system. 

The  matrix  equation  (3.41)  may  be  rewritten  as  a  set  of  2n  homo¬ 
geneous  equations, 


[♦  -  piK  =  o. 


(3.42) 
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Nontrivial  solutions  for  £  exist  if,  and  only  if,  the  inverse  of  the 
coefficient  matrix  [♦-pi]  does  not  exist  [36].  This  statement  leads 
to  the  characteristic  equation  of  $  given  by 

det  [♦  -  mU  =  0.  (3.43) 

This  expression  is  a  polynomial  equation  of  order  2n  in  the  scalar 
parameter  p.  Its  2n  roots  are  the  eigenvalues  of  $  and  may  be  real  or 
complex. 

When  the  elements  of  $  are  all  real,  it  is  easy  to  show  that 
complex  eigenvalues  occur  only  in  conjugate  pairs.  Taking  the  conju¬ 
gate  of  elements  of  both  sides  of  equation  (3.41)  results  in 

♦£  =  PC,  (3.44) 

where  ♦  is  unchanged  by  the  operation  since  its  elements  are  real,  and 
the  bar  above  the  remaining  elements  indicates  their  conjugate.  Writ¬ 
ing  the  characteristic  equation  of  ♦  from  equation  (3.44)  gives 

det  [♦  -  pi]  =  0.  (3.45) 

This  shows  that  the  conjugate  of  a  complex  eigenvalue  of  a  real  matrix 
is  also  an  eigenvalue  of  that  matrix.  As  a  result  the  complex  eigen¬ 
values  of  the  transition  matrix  must  occur  in  conjugate  pairs. 

Some  useful  relationships  between  any  matrix  and  its  eigenvalues 
may  be  obtained  by  comparing  the  characteristic  equation  of  the  matrix 
to  its  factored  form  in  terms  of  its  roots  (eigenvalues), 

=  (m1-p)(m2-m)  ...  (M2n-M)- 


det  [♦-pi] 


(3.46) 


When  p  is  set  equal  to  an  eigenvalue,  both  sides  of  equation  (3.46) 
are  zero  as  expected.  Letting  |J  =  0  gives  an  expression  for  the  deter¬ 
minant  of  the  transition  matrix, 

det  <t>  =  MaM2  •  •  •  M2n>  (3.47) 

i.e.,  the  determinant  of  a  matrix  equals  the  product  of  its  eigen¬ 
values.  Expanding  both  sides  of  equation  (3.46)  and  equating  like  co¬ 
efficients  of  p  gives  2n  relationships  including  (3.47).  From  the  co¬ 
efficient  of  p  to  the  2n-l  power,  it  is  determined  that  the  trace  of 
the  matrix  equals  the  sum  of  its  eigenvalues, 

2n  2n 

trace  $  =  1  <(>.  .  =  I  p.,  (3.48) 

i=l  i=l  1 

where  corresponds  to  the  diagonal  elements  of  4>. 

The  following  three  relationships  are  applicable  to  any  transition 
matrix  <t>  and  its  eigenvalues .  They  lead  to  an  important  property  of  the 
transition  matrix  of  a  Hamiltonian  system.  The  first  relationship  is 
that  the  eigenvalues  of  <J>  *  are  the  inverse  of  the  eigenvalues  of  <t>. 
This  is  shown  by  premultiplying  equations  (3.41)  by  4>  ,  dividing  the 

result  by  p,  and  rearranging  as  follows. 


The  second  relationship  is  that  eigenvalues  of  <t>  are  invariant  to 
a  similarity  transformation.  The  operation  C$C  where  C  is  a 
constant  matrix,  is  a  similarity  transformation  of  4>.  Premultiplying 
equation  (3.41)  by  C,  then  inserting  the  identity  expression,  C  1C,  in 
the  left-hand  side  of  the  result  verifies  this  relationship  as  follows 


(3.41) 

C4>£  =  Mce, 

OKfVo  =  m(c£). 

(3.50) 

The  third  relationship  is  that  eigenvalues  of  <l>  are  invar¬ 
iant  to  the  transpose.  This  may  be  shown  by  comparing  the  character- 

T 

istic  equation  (3.43)  for  ♦  with  that  for  its  transpose,  ♦  .  From  the 

construction  of  the  two  equations,  it  is  obvious,  by  interchanging 

row  and  column  operations  respectively  that  the  resultant  equations 

in  terms  of  p  are  identical.  Hence  the  eigenvalues  of  <l>  are  also  the 
T 

eigenvalues  of  . 

Applying  the  previous  result  (3.50)  to  this  result  and  choos¬ 
ing  the  matrix  C  to  be  the  fundamental  symplectic  matrix,  K,  gives 

KcV^Kt)  =  p(K£).  (3.52) 

When  4>  is  the  transition  matrix  of  a  Hamiltonian  system,  the  above 
equation  (3.52)  reduces,  using  its  symplectic  property  (3.32),  to 
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This  means  that  $  *  and  0  have  the  same  eigenvalues.  But  the  eigen¬ 
values  of  4>  1  were  also  shown  to  be  the  inverse  of  the  eigenvalues  of 
♦  .  This  implies  the  very  important  result  for  a  Hamiltonian  system, 
namely  that  the  eigenvalues  of  the  transition  matrix  4>  and  its  inverse 
must  occur  in  reciprocal  pairs. 

For  physically  realizable  systems  the  elements  of  the  transi' 
tion  matrix  are  real.  This  is  derivable  from  the  all  real  set  of  vari¬ 
ational  equations  determined  earlier  by  tacitly  assuming  a  physically 
realizable  system.  Combining  the  property  that  complex  eigenvalues  of 
real  valued  matrices  occur  in  conjugate  pairs  with  the  reciprocity 
property  of  eigenvalues  of  the  transition  matrix,  strong  restrictions 
on  the  complex  eigenvalues  of  a  transition  matrix  occur. 

A  single  pair  of  complex  eigenvalues  are  restricted  to  the 
unit  circle  in  the  complex  plane.  This  is  shown  by  first  applying  the 
reciprocity  property  to  a  pair  of  complex  eigenvalues,  where  j  =  V"l> 


'1 


=  a  +  jb, 


-  a  -  jb 

2  .  .2 

a  +  b 


(3.54) 

(3.55) 


In  order  that  and  also  be  conjugates,  the  denominator  in  (3.55) 
must  be  one.  Since  the  denominator  is  the  square  of  the  magnitude  of 
either  eigenvalue,  it  is  apparent  that  the  eigenvalues  lie  on  the  unit 
circle . 

For  more  complicated  systems,  the  complex  eigenvalues  may 
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occur  off  the  unit  circle  in  doubly  coupled  sets  of  four  eigenvalues. 
Each  set  is  composed  of  two  pairs  of  complex  conjugates,  which  may  also 
be  arranged  as  two  pairs  of  reciprocal  eigenvalues.  Many  other  config¬ 
urations  are  possible.  The  important  point  is  that  off  the  unit  circle 
complex  eigenvalues  are  restricted  to  occur  in  groups  of  at  least  four 
coupled  eigenvalues.  An  example  for  four  eigenvalues  would  be  the  re¬ 
lationships  (3.54)  and  (3.55)  plus  the  following  two, 


=  a  -  jb,  and 


(3.56) 


=  a_+Jb 

2  .  .2 

a  +  b 


(3.57) 


3 . 4  Properties  of  the  Monodromy  Matrix 

Many  of  the  properties  of  Hamiltonian  systems  with  periodic 
solutions  involve  the  monodromy  matrix.  This  special  constant  matrix 
maps  the  transition  matrix  forward  one  period  in  time.  After  deriving 
the  relationships  which  define  the  monodromy  matrix,  three  of  its 
important  properties  will  be  derived  to  conclude  this  chapter. 

Earlier  it  was  shown  by  construction  that  the  transition 
matrix  satisfies  the  variational  equation  (3.29).  Since  the  coeffi¬ 
cient  matrix,  A(t),  is  periodic  of  period  T,  then  4>(t+T,t0)  is  a 
solution  to  the  same  linear  differential  equations  as  4>(t,tQ).  From 
equation  (3.27),  the  two  solutions  must  be  related  by  a  constant 
matrix,  T,  as 

<Kt+T,t0)  =  O(t,t0)r.  (3.58) 


*■  '  j 
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Evaluating  equation  (3.58)  at  the  initial  time  defines  the  monodromy 
matrix  f, 


T  =  ♦(t0+T,t0).  (3.59) 

The  arguments  will  generally  be  dropped  for  the  monodromy  matrix, 
since  it  is  a  constant.  It  is  important  to  point  out  that  the  mono¬ 
dromy  matrix  varies  with  changes  in  the  initial  time  and  for  multiple 
periods.  These  relationships  will  be  established  in  the  two  following 
properties . 

The  first  property  of  the  monodromy  matrix  to  be  presented 
here  follows  directly  from  its  definition  in  equations  (3.58)  and 
(3.59).  The  relationship  of  the  monodromy  matrix  to  its  associated 
transition  matrix  is  determined  by  replacing  t,  in  equation  (3.58),  by 
t+T,  and  using  (3.58)  again  for  the  last  relation, 

*(t+2T,t0)  =  <Kt+T,t0)r  =  <Kt,t0)r2.  (3.60) 

Continuing  this  process  gives  the  general  formula, 

<J>(t+nT ,  tQ)  =  <Kt,tQ)rn,  (3.61) 

which  defines  the  functional  relationship  for  the  mapping  of  the 
transition  matrix  over  multiple  periods  by  the  monodromy  matrix. 

The  second  property  shows  that  even  though  the  monodromy 
matrix  is  dependent  on  the  initial  time,  t^,  its  eigenvalues  are  in¬ 
variant  to  such  changes.  Using  the  composition  rule,  the  relationship 
between  the  monodromy  matrix  initialized  at  some  time  t^  to  the  one  at 


-*  '  *■  v  ./ 
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some  other  time  t^  may  be  written  as 

♦(tj+T.tp  =  4>(t1+T,t2+T)*(t2+T,t2)<Kt2,t1), 

RT.tj)  =  *(t1+T,t2+T)r(T>t2)<Kt2,t1).  (3.62) 

Using  the  composition  rule  and  the  matrix  inverse  rule  for  the  first 
matrix  on  the  right  hand  side  of  (3.62)  gives 


♦  (tj+T.t^T)  =  «J>(t1+T,t0)[<Kt2+T,t0)] 


-1 


(3.63) 


From  the  definition  of  the  raonodromy  matrix  and  the  previous  rules, 
the  following  result  for  (3.63)  is  obtained 


♦  (tj+T.t^T) 


<Kt1,t0)r(T,t0)[<j»(t2,t0)r(T,t0)]  , 

*(t1,tQ)rr"V1(t2,t0), 

♦(tj,t2).  (3.64) 


Equation  (3.62)  can  now  be  expressed,  using  (3.64)  and  the  transition 
matrix  inverse  rule  as 


rcr.tj)  =  4."1(t2,t1)r(T,t2)<Kt2,t1).  (3.65) 

This  expression  shows  that  the  monodromy  matrix  for  a  Hamiltonian 
system  initialized  at  one  time  is  related  to  the  monodromy  matrix  of 
the  system  initialized  at  another  time  by  a  similarity  transformation. 
Thus  from  equation  (3.50),  the  eigenvalues  of  the  monodromy  matrix  are 


<  ■ 


..  .4  . 
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invariant  with  respect  to  its  initialization. 

The  third  property  shows  that  there  are  at  least  two  unity 
eigenvalues  of  the  monodromy  matrix  for  any  autonomous  Hamiltonian 
system  with  periodic  solutions.  Using  the  system  equations  expressed 
in  (3.9)  and  (3.10),  rewritten  below  for  convenience, 

y  =  Y(y),  (3.9) 

y(t+T)  =  y(t) ,  (3.10) 

the  variational  equations,  also  represented  by  equation  (3.25),  can  be 
expressed  as 

6y  =  6y  .  (3.66) 

Taking  the  time  derivative  of  equation  (3.9)  gives 


dy  . 

dt 


(3.67) 


which  shows  that  y(t)  is  a  solution  to  the  variational  equations  (3.66). 
From  equation  (3.27)  it  is  apparent  that  a  constant  vector  C  exists 
such  that 


y(t)  =  *(t,t0)C,  (3.68) 

where  <fr(t,tg)  is  the  transition  matrix  related  to  the  variational 
equations  (3.9).  In  other  words,  y(t)  must  be  a  linear  combination 
of  the  columns  of  the  transition  matrix.  At  the  initial  time  t  =  t^, 
the  transition  matrix  reduces  to  the  identity  matrix  and  the  constant 


/  L  '  >.  y  -y  ■  v 
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matrix  C  =  y(tQ).  Then  equation  (3.68)  may  be  rewritten 

y(t)  =  4>(t,t0)y(t0).  (3.69) 

After  one  period,  T,  equation  (3.69)  can  be  expressed  using  equations 
(3.9),  (3.10),  and  (3.68)  as 

y(t0)  =  4»(t0+T,t0)y(t0)  , 

y(tQ)  =  r(T,t0)y(t0).  (3.70) 

Rearranging  equation  (3.70)  in  the  form  of  equation  (3.40)  gives 

[r  -  I]y(t0)  =  0,  (3.71) 

which  shows  that  the  monodromy  matrix  has  one  unity  eigenvalue.  The 
eigenvector  corresponding  to  this  unity,  eigenvalue  is  the  tangent 
vector  to  the  periodic  path  of  the  solution.  It  is  defined  in  the  2n- 
phase  space  of  the  system  at  the  initial  time.  A  second  unity  eigen¬ 
value  of  the  monodromy  matrix  must  also  exist  due  to  the  previously 
reviewed  property  that  the  eigenvalues  of  a  Hamiltonian  system  must 
occur  in  reciprocal  pairs.  As  was  shown  in  the  first  section  of  this 
chapter  and  as  is  implied  by  the  two  unity  eigenvalues,  the  order  of 
the  system  is  reducible  by  at  least  two. 


CHAPTER  4 


NEW  THEORETICAL  RESULTS 

Contributions  to  the  general  theory  of  optimal  periodic 
control  developed  in  the  course  of  this  research  are  presented  in  this 
chapter.  A  principle  result  is  derived  in  each  of  the  following  four 
sections.  In  the  first  section,  a  new  algebraic  Riccati  equation  is 
derived,  which  determines  the  initial  conditions  for  generating  the 
periodic  solutions  of  the  Riccati  differential  equation  associated 
with  the  second  variation  of  the  performance  index.  An  extremely  use¬ 
ful  similarity  transformation  of  the  monodromy  matrix  is  developed 
from  this  result.  In  the  second  section,  a  large  class  of  extrema 
for  the  optimal  periodic  control  problem  is  shown  to  be  non-optimizing. 
The  symplectic  properties  of  the  monodromy  matrix  and  the  previous 
similarity  transformation  are  used  to  obtain  this  result.  A  weak 
sufficiency  condition  from  the  second  variation  is  developed  for  the 
free  period  case  in  the  third  section.  Finally,  a  periodic  regulator 
and  its  control  law  are  determined  by  investigating  optimal  paths 
near  a  local  periodic  optimum. 

4. 1  An  Algebraic  Matrix  Riccati  Equation 

A  new  algebraic  equation  is  derived  relating  the  Riccati 
variable  and  elements  of  the  monodromy  matrix.  The  variational  equa¬ 
tions,  which  determine  the  transition  matrix  for  the  state  variables 
and  Lagrange  multipliers,  are  obtained  from  the  Euler-Lagrange 
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equations,  (2.9)  through  (2.11),  by  taking  variations  about  an  extre¬ 
mal  solution.  Retaining  only  the  linear  terms  in  the  variations  gives 
the  following  2n  +  m  equations, 


fix  =  fx(t)6x  +  fu(t)fiu,  (4.1) 

6A  =  -H  (t)6x  -  H  (t)6u  -  fT(t)6A,  (4.2) 

xx  xu  x 

0  =  H  (t)6x  +  H  (t)6u  +  f^(t)6A>  (4.3) 

ux  uu  u 

where  the  partials  have  been  evaluated  along  the  extremal  path. 

Provided  Huu  is  nonsingular,  the  m  control  variations,  6u, 
may  be  written  from  (4.3)  in  terms  of  the  n  variations  in  the  state, 
fix,  and  the  n  variations  in  the  Lagrange  multipliers,  6A,  as 

6u  =  -(H  )_1 [H  fix  +  fT  6A)  .  (4.4) 

uu  1  ux  u 


Eliminating  the  control  variations  (4.4)  from  the  2n  equations  (4.1) 
and  (4.2)  determines  the  variational  equations  in  the  form  of  (3.25), 


fix 

.fiA. 


f  -f  (H  )  H 
x  u  uu  ux 


1  +  H  (H  )_1H 

XX  xu  uu  ux 


-f  (H  )_1fT 
u  uu  u 

■fix' 

T  -IT 

-f  +H  (H  )  f 

X  xu  uu  u 

,  (4.5) 


where  again  the  partial  derivatives  are  evaluated  along  the  extremal 
path.  Just  as  in  equation  (3.25),  the  elements  of  the  coefficient 
matrix  of  (4.5)  are  periodic  functions  of  time  with  the  same  period, 
T,  as  the  extremal  solutions. 
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The  transition  matrix  may  be  obtained  by  integrating  equation 
(3.29)  using  the  initial  condition,  •(tgjtg)  =  1  *  a&d  the  coefficient 
matrix  from  equation  (4.5).  This  results  in 


♦(t.to) 


(t,tQ) 

(t.tfl) 


* 


12 


* 


22 


(t,t0) 


(4.6) 


where  the  partitions  are  each  (n  x  n)  submatrices.  This  matrix  maps 
perturbations  in  the  state  and  Lagrange  multipliers  from  some  initial 
time,  tQ,  to  a  subsequent  time,  t,  as  follows. 


p5x(t)-| 

*U^fc,t0^  *12(t,V 

'6x(t0)' 

|_6A(t)  J 

*21  ^’*0^  *22^t,fc0^ 

6Mt0) 

An  expression  relating  a  perturbation  in  the  state  to  one 
in  the  Lagrange  multiplier  for  any  time  can  be  determined  by  manipulat¬ 
ing  equation  (4.7).  Assume  the  following  general  relationship. 


6A(t)  =  P(t)  6x(t), 


(4.8) 


where  P(t)  is  an  (n  x  n)  matrix  with  elements,  continuous  functions  of 
time.  Using  this  relationship  in  equation  (4.7)  to  eliminate  6\(tg) 
results  in, 


r6x(t)-i 

’*ii(t,t0)  +  ♦12(t»t0)P(t0)" 

L«Mt)J 

*21^t’t0^  +  *22^t,t0^t0^ 

-  *• 

which  for  convenience  may  be  abbreviated  as 


‘6x(t)‘ 

"X(t)‘ 

S 

.6A( t). 

L  ACtjJ  0 

6*(t0), 


(4.9) 


(4.10) 
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Provided  (X(t)]”1  exists,  the  variation  of  the  state  at  the 
initial  time,  (5x(tQ),  can  be  eliminated  in  equation  (4.10)  leaving  the 
result 

6A(t)  =  A(t)[X(t)]_16x(t).  (4.11) 

By  comparing  equations  (4.8)  and  (4.11),  the  following  expression  for 
P(t)  can  be  written, 

P(t)  =  A(t)[X(t)]'\  (4.12) 

which  in  expanded  form  is, 

P(t)  =  l<*>2l(t,t0)  +  *22(t,t0)P(t0)] 

•[4»u(t,t0)  +  ♦12(t,t0)P(t0)j"1.  (4.13) 

This  expression  for  P(t)  represents  the  solutions  to  the  Riccati-type 
differential  equation  (2.15),  rewritten  below,  which  was  derived  from 
the  second  variation, 

P(t)  =  I-fT+  H  (H  )_1fTJP(t)  -  P(t)[f  -  f  (H  )_1H  ] 

X  xu  uu  u  '  1  x  u  uu  UXJ 

+  P(t)[f  (H  )_1fT]P(t)  -  H  +  H  (H  )-1H  .  (2.15) 

U  UU  U  '  XX  XU  UU  UX 

The  details  of  showing  (4.13)  satisfies  (2.15)  are  presented  in 
appendix  A. 

The  transpose  of  P(t)  in  the  second  and  third  terms  of  equation 
(2.15)  has  been  dropped  since  P(t)  is  necessarily  symmetric.  This  is 
easily  shown  by  taking  the  transpose  of  both  sides  of  the  equation. 


-j" 
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T 

The  result  is  the  same  differential  equation  for  P  (t)  as  for  P(t), 
which  implies  that  P(t)  is  symmetric  if  it  was  initially  symmetric, 


P  (t)  =  P(t). 


(4.14) 


Now  assume,  that  a  periodic  solution  exists  to  the  Riccati 
type  equation  (2.15)  which  has  a  period,  T,  the  same  as  or  an  integral 
multiple  of  the  period  of  the  coefficients  of  P  in  the  equation.  An 
expression  for  the  initial  value  matrix  of  the  Riccati  variable,  P(t^), 
may  be  obtained  by  evaluating  equation  (4.13)  at  t  =  T  and  using  the 
periodicity  condition  P(t^+T)  =  P(tg).  This  provides 


P(t0)  I4>21(t0+T’V  +  <,*22^t0+T,t0^P^t0^ 


’(♦ll(tO+T»to)  +  <t,i2^t0+T,t0^P^t0^ 


(4.15) 


Using  the  relationship  (3.56),  between  the  monodromy  matrix 
and  the  transition  matrix,  in  the  above  expression  and  then  rearranging 
terms  results  in  the  following  algebraic  Riccati  equation 


P^t0^l'l2P^t0^  +  P<Vril  ”  r22P^t0^  “  r21  =  (4.16) 

Solutions  to  this  equation  provide  the  initial  conditions,  P(tg),  that 
generate  all  of  the  periodic  solutions  to  the  Riccati  differential 
equation  (2.15). 

Since  the  initial  time,  t^,  is  arbitrary,  the  algebraic 
expression  (4.16)  is  valid  for  any  time.  The  time  relationships 
for  the  coefficients  of  P  are  determined  using  the  similarity  trans¬ 
formation  (3.63)  when  the  monodromy  matrix  is  known  for  some  initial 
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tine,  tj,  and  period  T.  The  coefficients  are  then  the  appropriate 
blocks  of  the  resultant,  time-varying  monodrony  matrix, 

T(T,t)  =  *(t,t1)r(T,t1)*“1(t,t1),  (4.17) 

where  the  time  tg  has  been  replaced  by  t  to  emphasize  the  variability 
in  time.  In  this  manner,  the  Riccati  variable  may  be  determined  for 
any  time  using  a  given  monodromy  matrix  and  the  transition  matrix. 

Another  useful  result  of  equation  (4.16)  is  derived  from  a 
similarity  transformation  of  the  monodromy  matrix.  An  off-diagonal 
block  of  the  transformed  matrix,  partitioned  into  four  (n  x  n)  blocks, 
is  reduced  identically  to  zero.  The  relationships  and  properties 
derived  from  this  simple  operation  are  so  important  to  the  new  develop¬ 
ments  in  this  chapter  that  the  transformation  is  presented  here  as  a 
theorem. 


THEOREM  4.1:  If  a  periodic  solution  to  the  matrix  Riccati 
differential  equation  (2.15)  exists  with  period  T  (the  same  or  an 
integral  multiple  of  the  period  of  the  coefficients  of  P  in  the  equa¬ 
tion),  then  the  monodromy  matrix,  given  by 


!~(T,t0)  = 


F11  r!2 
r21  r22 


(4.18) 


is  similar  to  the  transformed  matrix 


-1 


.  [rntr12 


P(t0> 


■  12 

r22-P(Vri2 


Lrcr,t0)L 


o 


(4.19) 
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where  the  similarity  transformation  matrix,  L,  is  defined  by 


L  = 


'  I 


(4.20) 


and  P(tg)  is  the  initial  value  matrix  that  satisfied  the  algebraic 
Riccati  equation  (4.16). 


PROOF:  Perform  the  matrix  operations  of  the  similarity 
transformation  indicated  by  the  left-hand  side  of  (4.19)  where  the 
inverse  of  L  is 


p(t0> 

This  results  in  a  matrix  with  the  three  non-zero  blocks  of  (4.19)  and 
a  term  which  is  the  algebraic  Riccati  equation  (4.16)  in  the  remaining 
block.  Due  to  the  choice  of  P(t^) ,  this  block  is  identically  zero  and 
the  transformed  matrix  is  the  desired  one. 


O' 

I. 


(4.21) 


4.2  Eigenvalues  and  Sufficient  Conditions 

It  was  shown  in  section  3-2  that  the  monodromy  matrix  of  the 
Hamiltonian  system  defined  by  the  periodic  control  problem,  (1.1) 
through  (1.3),  is  symplectic.  The  resulting  properties  of  the  mono¬ 
dromy  matrix  and  those  associated  with  the  transformed  matrix  of 
theorem  4.1  from  the  previous  section,  restrict  the  eigenvalues  of  the 
monodromy  matrix  for  optimal  solutions.  It  is  shown  that,  except  for 
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possibly  isolated  limit  points,  local  optimal  solutions  are  all  un¬ 
stable.  The  sufficient  conditions  [14]  for  a  local  minimum,  presented 
in  section  2.2,  are  thus  modified  by  the  developments  in  this  section. 

Since  the  monodromy  matrix  is  symplectic,  it  must  satisfy  the 
relationship  (3.14)  resulting  in  the  following  equation, 

rxrT  =  K,  (4.22) 


where  K  is  the  fundamental  symplectic  matrix  (3.7).  Writing  this  rela¬ 
tionship  in  block  form,  consistent  with  the  partitioning  of  equation 
(4.18),  and  using  equations  (3.17)  through  (3.19)  gives 


r  rT 

11  12 


T 

r  r 

22  21 


r  rT 

'  11  22 


T 

r  r 

1 12  11  * 


=  r  r 

'  21  22  » 


r  r  =  I 
1221 


(4.23) 

(4.24) 

(4.25) 


The  symplectic  property  of  the  monodromy  matrix  is  preserved 
through  the  similarity  transformation  in  the  theorem  of  the  preceding 
section.  This  may  be  demonstrated  by  showing  that  the  transformed 
matrix  LIT  1  satisfies  the  following  relationship  derived  from  the 
symplectic  property  (3.14), 

(LI~L~ 1  )K(Ln" 1  )T  =  K.  (4.26) 

-1  -T 

Premultiplying  both  sides  of  (4.26)  by  L  and  postmultiplying  by  L 
gives 


r(L’1KL'T)rT 


(l_1kl‘t). 


(4.27) 
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If  the  matrix  L  is  symplectic,  then  from  (3.14)  the  following  equiv¬ 
alent  relationships  hold, 

LKLT  =  K  and  K  =  L_1KL~T,  (4.28) 

and  equation  (4.27)  reduces  to  the  identity  given  by  (4.22).  Perform¬ 
ing  the  matrix  operations  indicated  by  either  of  equations  (4.28) 
verifies  that  both  matrices  L  and  Ln  *  are  symplectic. 

Expressing  LrL  1  in  the  block  form  (4.19)  and  expanding 
equation  (4.26)  to  obtain  relationships  between  combinations  of  ele¬ 
ments  of  the  resulting  symplectic  matrix,  as  in  (3.17)  through  (3.19), 
provides  only  one  new,  but  very  important  equation, 

[rn  +  r!2P||r22  -  Pri2)T  =  ’■  01 

I r 22  -  Pri2]T  =  (ril  +  ri2P]1  •  (*’29) 

The  significance  of  this  new  equation  is  that  it  strongly 
restricts  the  eigenvalues  of  the  monodromy  matrix  corresponding  to 
real-valued  Riccati  variable  elements.  Note  that  the  eigenvalues  of 
the  two  matrices,  [F^  +  r^PJ  and  ~  PF^l  from  equations  (4.19) 

and  (4.29),  are  also  the  eigenvalues  of  the  monodromy  matrix,  F,  since 
eigenvalues  are  invariant  through  a  similarity  transformation,  as  estab 
lished  previously  by  (3.50).  The  transformation  also  assigns  the  eigen 
values  of  the  monodromy  matrix  to  the  two  submatrices  on  the  diagonal 
of  the  transformed  matrix  such  that  the  eigenvalues  of  one  are  the 
reciprocal  of  those  of  the  other  submatrix. 
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By  considering  small  perturbations  about  the  extremal 
solution,  the  restrictions  on  the  eigenvalues  due  to  equation  (4.29) 
become  apparent.  The  monodromy  matrix  maps  perturbations  over  one 
complete  cycle  or  period.  The  elements  of  this  matrix  (a  special  case 
of  the  transition  matrix)  are  always  real-valued,  provided  the  applica¬ 
tion  is  physically  realizable.  Also,  the  elements  of  the  associated 
Riccati  matrix  must  be  real-valued  for  the  extremal  solutions  to  be 
locally  optimizing.  As  a  result  each  element  of  the  two  matrices  of 
(4.29)  must  be  real-valued,  requiring,  from  (3.44)  and  (3.45),  that 
the  complex  eigenvalues  of  each  matrix  must  occur  in  conjugate  pairs. 
The  reciprocal  relationship  (4.29)  between  the  two  matrices  requires 
that  complex  eigenvalues  of  the  monodromy  matrix,  on  the  unit  circle, 
must  occur  in  even  multiple  conjugate  pairs,  i.e.,  double  conjugate 
pairs,  quadruple  conjugate  pairs,  and  so  forth.  Generally,  this  spe¬ 
cial  circumstance  results  at  a  limit  point  where  coupled  eigenvalues 
coalesce  upon  entering  or  exiting  the  unit  circle. 

Periodic  solutions  are  commonly  classified  as  stable  when  all 
eigenvalues  of  the  monodromy  matrix  lie  on  the  unit  circle  and  unstable 
when  at  least  one  pair  lies  off  the  unit  circle.  Consequently,  locally 
optimizing  periodic  solutions  are  unstable  with  the  possible  exception 
of  isolated  limit  points.  This  statement  further  restricts  the  suffi¬ 
cient  conditions  [14]  for  a  local  optimum  presented  in  section  2.2. 

A  deficiency  in  the  original  statement  of  these  conditions  is 
now  examined.  Recall  from  section  3.4  that  the  monodromy  matrix  for  a 
periodic  Hamiltionian  system  always  has  at  least  two  unity  eigenvalues. 


Due  to  the  relationship  (4.29),  at  least  one  unity  eigenvalue  will 
occur  in  each  of  the  submatrices  on  the  diagonal  of  the  transformed 
matrix  (4.19).  Since  Z  of  equation  (2.17)  is  identically  equal  to 
(rn  +  F 22^^  *  one  t*le  submatrices  of  (4.19),  the  condition  (5)  of 
the  sufficiency  conditions  is  never  satisfied. 

From  equation  (3.71)  the  eigenvector  corresponding  to  one  of 
the  unity  eigenvalues,  is  tangent  to  the  reference  (extremal)  solution 
in  phase  space.  As  a  result,  small  perturbations  along  the  tangent 
path  do  not  change  the  amplitude  or  period  of  the  perturbed  solution. 
The  only  apparent  effect  when  comparing  the  perturbed  solution  to  the 
unperturbed  solution  is  a  time  or  phase  shift.  If,  as  in  the  general 
problem  (1.1)  through  (1.3),  the  only  condition  is  the  periodicity 
constraint,  the  optimal  periodic  solution  is  invariant  to  a  perturba¬ 
tion  along  its  tangent. 

4.3  Free  Period  Second  Variation  Condition 

The  weak  sufficient  conditions  for  optimality  derived  by 
Bittanti,  Locatelli,  and  Maffezzoni  [14]  for  periodic  control  problems 
with  fixed  period  are  extended  in  this  section  to  the  class  of  problems 
for  which  the  period  is  unrestricted.  The  modified  conditions,  two  of 
which  are  also  applicable  to  the  fixed  period  case,  are  developed 
from  the  second  variation  applying  results  obtained  from  the  Hamilton- 
Jacobi  theory. 

The  second  variation  may  be  expressed  in  terms  of  operations 
on  the  first  variation,  given  by  equation  (2.8)  as  follows, 
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It=T 

}  +  d{|(H-J)dT} 

+  d{  |\T1(H  ♦  XT)6x  +  H  6u  ♦  6\T(hT  -  x)]dt} .  (4.30) 

x  u  A 


The  first  variation  conditions,  (2.9)  through  (2.14),  must  be  satisfied 
as  well  as  the  following  variational  constraints, 


6x  =  f  6x  +  f  6u,  (4.31) 

xu* 

dx(T)  =  dx(0) ,  (4.32) 


which  are  the  variations  of  the  system  constraints,  (1.2)  and  (1.3), 
to  the  first  order. 

The  expression  for  the  second  variation  (4.30)  may  be  reduced 
in  the  following  manner.  Using  equations  (1.3),  (2.13),  and  (4.32)  in 
the  first  term  of  the  expression  gives, 


d{dvT4*  +  (vT-|;XT)dx 


t=T 

} 

t=0 


TIT  1  T 
(2dv1+i  X  dT~  dX  )dx 


t=T 

t=0, 


1  ..T. 
jdX  dx 


t=T 


t=0 


(4.33) 


Using  equations  (2.11),  (2.14),  and  dJ  =  0  in  the  second  term  of  the 
expression  (4.30)  gives 


d{|(H-J)dT}  =  -^(H-J)dT2  +  £  (dH-dJ)dT, 

I  t=T 


|(Hxdx  +  dXTH*) 


dT  . 


(4.34) 
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i  the  last  term  of  expression  (4.30),  using  conditions  (2.9)  through 

•T 

2.11),  equation  (4.31),  and  integrating  by  parts  the  6\  term  gives 


d{|  1  [ELjdt}  =  -- 
JO  T 


=  >i 


T  fT 

[EL]dt  +  i[EL]dt  +  H  6[EL]dt, 

0  t=0  Jo 


=  |(6AT6x) 


t=0  T  3  0 


[6xT6uT]  |Hxx  Hxulr6xl  dt,  (4.35) 


H  H  I  6u 
ux  uu  *- 


where  (EL]  =  ( (H  +\T)6x  +  H  6u  +  6A^(h7-x)]. 

X  U  A 

Combining  these  results,  equations  (4.33)  through  (4.35), 
for  the  second  variation  provides  the  following  expression, 

t=T  I t=T 

d2J  =  i((6AT6x  -  d\Tdx)  +  (H  dx  +  d\TH^)l  dT] 

1  t=0  x  A  1 


W6uT]  H*x  Hxu  PXl 


H  H  I 6u| 
ux  uu  I1-  J 


(4.36) 


Replacing  the  6A  6x  term  using  the  relationship  (1.4)  and  conditions 
(2.9)  through  (2.11)  gives  a  more  useful  form  of  the  second  variation, 


d2J  =  |(dxT(T)dT] 


H*(T)  irdxd)-] 

-Hx(T)f(T)JL  dT  J 

H  H  |  f  6x  1 

XX  XU  dt, 

H  H  I  6uJ 
ux  uu  J  u 


Hx(T)  -  Hx(T)f(T) 


T(6xT6uT] 


(4.37) 


subject  to  the  variational  constraints  given  by  (4.31)  and  (4.32) 
As  in  the  accessory  minimum  problem  [151,  showing  that 
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d  J  >  0  for  any  6u(t)  t  0  establishes  the  minimality  of  the  extremal 
solution.  This  may  be  accomplished  by  adding  to  equation  (4.37)  the 
identically  zero  quantity, 


S6x^P(f  fix  +  f  6u  -  6x)dt  =  0, 
0 


(4.38) 


where  the  partial  derivatives  of  f  and  H  are  evaluated  along  the  ex¬ 
tremal  trajectory,  and  P  is  the  Riccati  variable  determined  by  the 
matrix  differential  equation  (4.14).  The  integral  of  the  resulting 

relationship  for  the  second  variation  can  be  written  as  a  perfect 

« 

square  after  integrating  the  fix  term  by  parts  and  using  equation 
(4.14)  to  appropriately  express  the  coefficients  of  the  variations. 
This  gives  for  the  second  variation 


d2J  =  i  {2Hx(T)dx(T)dT  -  Hx(T)f(T)dT2-  (fixT(t)P(t)fix(t)) 


|t=T 


t=0 
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[fixT  fiuT] 


(H  +Pf  )  (H  )"1fTP+n  )  (H  +Pf  )' 
XU  u  uu  u  ux  ux  u 
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uu 


(4.39) 


which  is  commonly  expressed  in  the  norraed  form  as, 


d2J  =  |  {2Hx(T)dx(T)dT  -Hx(T)f (T)dT2-(6xT(t)P(t)6x(t)) 


.i(T 

T)o 


|(H  )_1(fT  P  +  H  ) fix  +  5i 

1  uu  u  ux 


“1|h  jt- 

1  UU 


t=T 

} 

t=0 

(4.40) 


The  end  point  expression  in  (4.40)  will  be  shown  to  vanish.  First  the 
following  relationship  must  be  developed. 


MO)  =  \(0)  -  P(0)  x  (0)  =  0. 


(4.41) 


From  equation  (3.71),  the  relationship  between  the  monodromy 


matrix,  a  unity  eigenvalue,  and  the  associated  eigenvector  may  be 
expressed  as 

Ty(0)  =  y(0) ,  (4.42 

where  the  eigenvector  y(0)  is  the  phase  space  velocity  vector  given 
by 

r*(°n 

y(0)  =1  ,  (4.43 

LMO)J 

Performing  a  transformation  of  variables  in  (4.42)  using  the 
similarity  transformation  matrix,  L,  defined  by  (4.20)  provides 
the  following  relationship, 

(LrL_1]Ly(0)  =  Ly(0)  .  (4.44 

By  substituting  in  (4.44)  the  equivalent  expressions  given  by 
(4.19),  (4.41),  and  (4.43),  the  following  matrix  equation  is  obtained, 


Assume  first  that  \  ^  0  and  satisfies  the  bottom  equation 

of  (4.45),  which  may  be  expressed  as 

» 

IX22  '  PA12  '  11  A  =  0  *  (4.46] 

_T 

Premultiplying  the  upper  equation  of  (4,45)  by  A  gives  the  result, 
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*  •  • 

XT(ru  +  r12P  -  ni  +  XTr12x  =  o.  (4.47) 

The  first  term  in  this  expression  is  reduced  to  zero  by  the  following 
operations.  First,  multiply  equation  (4.46)  by  (I~22  -  Pl"12)  \ 
which  gives 

II  -  (r22  -  pr12)_1]x  =  0.  (4.48) 

Using  the  symplectic  property  derived  in  (4.29)  for  the  transformed 
monodromy  matrix  provides 

[I  "  (ril  +  ri2P)T]X  =  °*  (4-49) 

which  after  transposing  becomes 


xVn  +  r12P  -  i)  =  o. 


(4.50) 


Using  this  expression  reduced  equation  (4.47)  to 


xTr  12x  =  0. 


(4.51) 


This  equation  may  be  written  in  terms  of  the  symmetric  part  of  f"^, 


xt  ( r>2/  r'>y>. 


(4.52) 


since  the  result  is  a  scalar.  However,  in  general 

r-T 


* 


X  *  0 


(4.53) 


for  \  f  0,  obtained  from  equation  (4.46),  which  contradicts  the 


original  supposition;  thus 


A  =  0,  (4.54) 

justifying  the  relationship  (4.41). 

Using  the  relationship  for  the  variations  (1.4),  the  rela¬ 
tionships  (2.9)  through  (2.10),  and  the  result  (4.41);  the 
end  point  expression  in  (4.40)  is  identically  zero,  as  shown  below, 

2H«dXdT  -  Hx£°dT2|t=T  *  6*TpHt=T  -  S*TPC*|t=0. 

=  ~(xTPdx  +  dxTPx  -  xTPx  dT)dT 

t=T 

=  2H°dxdT  -  H°f°dT2.  (4.55) 

This  leaves  for  the  second  variation  the  squared  expression, 

“2j  *  f J]K„)''(fuTp  +  o*  * «» |  [o dt-  <<"56) 

uu 

The  neighboring  optimal  control  law  that  causes  the  second  variation 

(4.56)  to  vanish  is 

6u°  =  -(H°  )_1(f0  P  +  H°  )  6x.  (4.57) 

uu  u  ux  v  * 

Any  other  control  variation  produces  a  greater  cost  since  from  equation 

(4.56) ,  d2J  £  0. 

In  order  to  show  that  the  extremal  solution  is  locally  mini¬ 
mizing,  it  must  be  shown  that  there  is  no  neighboring  periodic  solution 
for  which  the  generating  variations  satisfy  the  boundary  condition 
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(4.32)  and  use  the  control  law  (4.57).  The  boundary  condition  may  be 
rewritten  using  equation  (1.4)  in  (4.32)  and  recognizing  that  f°  is 
,  periodic  with  period  T, 

6x(0)  =  6x(T)  +  £°(0)dT.  (4.58) 

Using  the  neighboring  optimal  control  law  (4.57)  in  the  dynamic  con¬ 
straint  equation  (4.31)  provides  the  following  linear,  homogeneous 
differential  equation  for  propagating  the  variation  in  the  state, 

6k  =  (f°  -  f°(H°  )_1(f°TP  +  H°  ))  6x.  (4.59) 

X  u  uu  u  ux 

This  equation,  (4.59),  is  the  same  as  the  first  of  the  set 
of  equations  (4.5)  with  6A  replaced  by  (4.8).  The  corresponding 
transition  matrix  from  (4.9)  may  be  expressed  as 

4>z  =  4*^+  4>12P,  (4.60) 

where  Z  is  the  coefficient  matrix  of  (4.59),  also  defined  in  (2.17). 

It  was  shown  in  section  4.2  that  the  transition  matrix  4>z,  evaluated 
over  one  period  and  expressed  in  terms  of  elements  of  the  monodromy 
matrix,  is 

4>Z(T,0)  =  rn  +  ri2P(0),  (4.67) 

and  has  at  least  one  unity  eigenvalue.  Relative  to  this  eigenvalue  is 
the  eigenvector  f°(0),  a  vector  tangent  to  the  optimal  periodic  path 
at  the  initial  time.  This  is  the  final  relationship  developed  at  the 
end  of  chapter  three. 
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The  class  of  all  possible  neighboring  solutions  is  considered 
by  examining  the  initial  variation,  6x(0),  that  generates  them.  A 
variation  parallel  to  the  eigenvector  f°(0)  and  an  other  one  orthogonal 
to  it  are  examined.  Of  concern  is  whether  or  not  a  set  of  variations 
can  be  found  for  which  the  periodicity  constraint  (4.58),  is  satisfied 
while  the  control  law,  (4.57),  through  the  transition  matrix  (4.61)  is 
applied.  Any  initial  variation,  6x(0),  can  be  represented  by  a  combi¬ 
nation  of  these  two  variations. 

The  relationship  between  the  initial  variation  and  the  vari¬ 
ation  corresponding  to  the  extremal  (reference)  trajectory  after  one 
period  is  given  by 

6x(T)  =  4>z(T,0)6x(0).  (4.62) 

When  the  initial  variation,  6x(0),  is  parallel  to  f°(0),  the  relation¬ 
ship  (4.62)  reduces  to 

6x(T)  =  6x(0).  (4.63) 

Applying  this  result  to  the  constraint  (4.58)  implies  that  the  period 
is  unchanged  since  the  initial  time  is  arbitrary,  and  a  non-zero  f°(0) 
can  always  be  chosen.  The  variation  at  the  initial  time  represents  a 
change  from  one  point  on  the  extremal  path  to  another  point  on  the 
extremal  path.  Consequently,  the  component  of  a  variation  parallel  to 
f°(0)  produces  no  change  in  the  optimal  trajectory. 

When  the  initial  variation  is  orthogonal  to  f°(0),  it  is 
found,  by  applying  the  relationship  (3.42)  to  each  component  of  6x(0) 
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relative  to  the  eigenvectors  (basis  vectors)  of  that  fix(T)  is  also 
orthogonal  to  f°(0).  Provided  the  eigenvalues,  excluding  the  unity 
eigenvalue  corresponding  to  f°(0),  do  not  lie  on  the  unit  circle,  the 
following  relationship  is  true  for  any  variation  orthogonal  to  f°(0), 

6x(T)  t  6x(0).  (4.64) 

However,  because  neither  variation  has  a  component  in  the  direction  of 
f°(0),  there  is  a  contradiction  between  (4.58)  and  (4.64).  Therefore, 
with  the  conditions  on  the  eigenvalues  of  ^  stated  above,  no  neighbor¬ 
ing  solution  can  be  found  which  satisfies  (4.58)  and  uses  the  control 
law  (4.57).  As  a  result,  the  extremal  solution  is  a  local  minimum. 

The  sufficient  conditions  of  Bittanti,  Locatelli,  and 
Maffezzoni  [14]  for  local  optima  of  periodic  control  problems  as  mod¬ 
ified  and  developed  in  this  section  are  summarized  in  the  following 
theorem  for  either  specified  or  unrestricted  period: 

THEOREM  4.2.  For  the  periodic  control  problem  defined  by 
equations  (1.1),  (1.2),  and  (1.3),  the  control  u°(*)»  a  piecewise- 
continuous,  periodic  function  of  the  input  space,  is  a  local  minimum 
if  the  following  conditions  are  satisfied: 

(1)  The  Euler-Lagrange  equations,  (2.9)  through  (2.11), 

and  the  transversality  condition  (2.13)  are  satisfied, 
where  u°(*)  is  a  solution  to  (2.11)  for  0  £  t  £  T  and 
x°(«)  and  A°(*)  are  the  corresponding  solutions  to  the 
two-point  boundary  value  problem  which  results  from 
introducing  u°(*)  into  (2.9)  and  (2.10).  When  the 
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period  is  unrestricted,  the  additional  transversality 
condition  (2.14)  must  also  be  satisfied; 

(2)  The  strong  form  of  the  Legendre  condition  is  satisfied, 

H°  >  0,  for  0  S  t  £  T; 
uu 

(3)  A  bounded  symmetric  solution  to  the  Riccati  equation 
(2.15),  exists  for  0  i  t  S  T  subject  to  the  periodic 
boundary  condition,  (2.16);  and 

(4)  Except  for  the  unity  eigenvalue  associated  with  the 
eigenvector  f°(0),  no  eigenvalue  of  ^(TjO)  lies  on 
the  unit  circle. 

4.4  A  Periodic  Regulator 

In  the  previous  sections,  an  open-loop  solution  to  the  opti¬ 
mal  periodic  control  problem  (1.1),  (1.2),  and  (1.3)  was  developed. 
Examining  the  restrictions  on  the  eigenvalues  of  the  monodromy  matrix, 
specified  in  section  4.2,  it  is  apparent  that  except  for  very  special 
circumstances,  the  optimal  open-loop  solution  is  unstable.  In  this 
section  a  control  strategy,  analogous  to  the  static  regulator,  is 
presented  which  minimizes  the  cost  of  holding  a  perturbed  periodic 
system  near  its  optimal  open-loop  trajectory. 

Perturbations  of  a  local  optimal  trajectory  may  be  expressed 
in  terms  of  the  variational  equations  (4.5).  Using  the  transformation 
matrix  (4.21)  of  Theorem  4.1,  a  new  set  of  variations  may  be  defined  as 


"  6x(t)  -| 

-  i  on 

-6x(t)- 

.  6A*(t)J 

.-p(t)  ij 

.6Mt). 

(4.65) 
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where  P(t)  is  a  solution  of  the  Riccati  equation  (4.14).  In  terms  of 
these  transformed  variables,  the  variational  equations  become 


6x(t) 

"  Z  -f°(H  ®)-1f®T" 
u  uu  u 

6x(t) 

6A*(t)_ 

T 

0  -Z 

6A*(t) 

(4.66) 


where  Z  is  defined  by  equation  (2.17)  and  the  zero  element  is 
ity  expression  of  the  Riccati  equation  (4.14).  Note  that  the 
formed  variation,  6A*(t),  is  identically  zero  due  to  equation 
the  relationship  between  variations  of  the  state  and  co-state 
multipliers).  Therefore  equation  (4.66)  reduces  to 


an  ident- 
trans- 
(4.8)  and 
(Lagrange 


6x(t)  =  Z  6x(t), 


(4.67) 


where  the  boundary  condition  6x(0)  is  arbitrary.  This  equation  propa¬ 
gates  a  perturbation  to  an  optimal  trajectory  using  the  control  law 
(4.57),  which  was  incorporated  into  the  variational  equations  through 
the  transformation  (4.65). 

For  the  linear  regulator  problem,  as  in  Bryson  and  Ho  [15], 
the  coefficients  in  the  variational  equations  (4.5)  are  constants. 

This  is  the  case  for  the  static  equilibrium  solutions  associated  with 
the  periodic  control  problem.  In  addition,  the  Riccati  variable  in 
(4.65)  is  a  steady  state  solution  of  the  Riccati  equation  (4.14). 
Therefore  Z  in  (4.55)  is  also  a  constant.  Provided  the  system  is  con¬ 
trollable,  then,  as  in  Brockett  [36],  a  solution  to  P  =  0  exists  such 
that  the  real  parts  of  the  eigenvalues  of  Z  are  negative  and  6x(t)  ■*  0 
as  t  ■>  ®.  Using  this  P  in  the  control  law,  equation  (4.59),  provides 


the  relationship  for  the  static  regulator  of  a  linear  system, 


6u° 


■(H°  )-1foTP6x. 
uu  u 


(4.68) 


Consider  now  a  regulator  for  an  optimal  periodic  trajectory 
(orbit).  As  indicated  in  section  4.1,  the  coefficients  of  the  vari¬ 
ational  equation  (4.5)  are  periodic.  The  relationship,  (4.7),  for  the 
transition  matrix  shows  how  a  perturbation  will  propagate  away  from 
an  optimal  orbit  in  the  open  loop  control  no  feedback  case.  Provided 
a  periodic  solution  to  the  Riccati  equation  (4.14)  exists,  Theorem  4.1 
shows  that  in  terms  of  the  transformed  variables  (4.65)  the  propagation 
equations  (4.7),  evaluated  over  one  period,  T,  becomes 


'fix(T) 

.5\*(T) 


ril  +  ri2P  ri2  1 

’6x(0) 

0  "22  -  Pr!2 

J 

6A*(0)_ 

(4.69) 


The  initial  variation,  6x(0),  for  (4.57)  is  arbitrary,  and  as  before 
the  transformed  variation,  6A*(t),  is  identically  zero.  Therefore,  the 
discrete  equation  (4.69)  may  be  reduced  for  multiple  stages  to 


6x(nT)  =  [Tn  +  ri2Pj  6x((n-l)T) , 

=  <t>z  6x((n-l)T), 

=  n  6x  (0),  (4.70) 

since  the  matrix  [1~jj  +  r^PJ  is  a  constant  from  stage  to  stage,  its 
components  having  been  evaluated  at  the  time  of  the  perturbation.  As 
in  equation  (4.67)  the  optimal  control  law  (4.57)  is  incorporated  in 
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the  discrete  equation  (4.70)  which  propagates  perturbations  to  the 
optimal  orbit. 

In  the  closed  loop  (4.58)  the  optimal  orbit  forms  a  natural 
limit  cycle  for  state  perturbations  provided  that  the  eigenvalues 
of  <J>£  satisfy  Theorem  4.2,  condition  (5).  The  components 
and  ri2  of  <t>2  in  (4.70)  are  both  elements  of  the  monodromy  matrix 
of  the  optimal  orbit  and  are  uniquely  determined  by  the  initial  time 
associated  with  the  arbitrary  state  perturbation,  6x(0).  The  remain¬ 
ing  component,  P,  is  one  of  the  solutions  to  the  new  algebraic 
Riccati  equation  (4.17).  Each  solution  represents  an  initial  condition 
which  produces  a  periodic  solution  to  the  Riccati  differential 
equation  (4.14).  As  a  result  there  is  a  one  to  one  correspondence 
between  solutions,  P,  to  the  algebraic  Riccati  equation  and  permis¬ 
sible  combinations  of  the  eigenvalues  of  the  monodromy  matrix, 
which  comprise  the  eigenvalues  of  Therefore,  the  periodic  regu- 

ator  is  the  control  law  (4.57)  determined  by  the  particular  solution, 

P,  to  (4.17)  for  which  the  eigenvalues  of  satisfy  Theorem  4.2, 
condition  (5).  In  some  respects  this  is  a  generalization  of  the 
linear  regulator  problem. 

A  concept  for  implementing  the  periodic  regulator  may  be 
developed  which  exploits  two  properties  of  the  regulator.  First,  the 
gains  for  the  regulator  are  dependent  on  parameters  of  the  reference 
(optimal)  orbit  only.  The  time  of  the  perturbation  may  be  identified 
by  a  point  on  the  orbit.  Therefore,  the  gains  may  be  computed  apriori 
relative  to  states  on  the  reference  orbit  and  stored  in  a  look-up 
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table.  Second,  the  closed  loop  periodic  system  is  invariant  to  pertur¬ 
bations  tangent  to  the  reference  orbit.  As  a  result,  from  an  orthog¬ 
onal  projection  of  the  perturbed  state  onto  the  reference  orbit,  both 
the  new  perturbation  and  the  corresponding  control  gains  can  be  deter¬ 
mined  for  any  desired  interval  along  the  perturbed  path. 


CHAPTER  5 


AN  OPTIMAL  PERIODIC  CONTROL  PROBLEM 

A  particular  optimal  periodic  control  problem  is  formulated 
in  this  chapter  and  some  preceding  results  are  simplified  by  exploit¬ 
ing  the  order  and  symmetry  of  the  problem.  The  dynamics  are  chosen  to 
be  simple,  consisting  of  a  double  integration  of  the  scalar  control 
variable.  The  performance  index,  which  is  not  convex  in  the  state 
variables,  is  constructed  to  allow  optimal  periodic  solutions  for  a 
range  of  control  variable  weighting.  The  principal  objective  of 
this  chapter  is  to  develop  an  example  problem  for  which  a  local 
optimal  periodic  solution  may  be  obtained  that  illustrates  some  of 
the  general  characteristics  of  the  class  of  optimal  periodic  control 
problems.  Recent  efforts  are  briefly  reviewed  in  the  first  section  to 
provide  the  background  from  which  the  present  research  effort  is 
motivated.  In  the  next  section  the  example  problem  is  formulated,  the 
first  order  necessary  conditions  are  examined,  and  the  frequency  test 
is  applied.  Some  simplifications  of  results  presented  in  previous 
chapters  are  developed  in  section  three  for  fourth  order,  symmetric 
systems  such  as  the  example  problem.  Finally,  in  the  last  section, 
an  explicit  solution  to  the  Riccati  equation  is  derived  for  any  sym¬ 
metric,  fourth  order,  periodic  system. 

5 . 1  Background 

An  interesting  application  of  optimal  periodic  control  is 


developed  in  a  lengthy  controversy  [5,  6,  7,  8]  over  a  particular 
aerospace  problem.  The  contention  pertains  to  the  optimality  with 
respect  to  fuel  consumption  of  the  steady-state  cruise  of  an  aircraft. 
The  results  obtained  by  Speyer  (5]  stimulated  additional  work  [9], 
attempting  to  find  an  optimal  periodic  solution  for  the  cruise  segment 
of  an  aircraft  flight  path.  Some  of  the  important  aspects  of  this 
aerospace  problem  are  now  briefly  reviewed  to  establish  the  background 
and  rationale  for  constructing  a  new  and  simpler  optimal  periodic  con¬ 
trol  problem.  The  fabricated  problem  is  thoroughly  examined  in  the 
remainder  of  this  work. 

Steady-state  cruise  is  established  for  an  aircraft  by  mini¬ 
mizing  its  range  factor,  defined  as  fuel  rate  per  range  rate.  All 
control  variables  are  held  constant  when  operating  at  this  static 
cruise  point.  A  frequently  used  dynamic  model  for  evaluating  aircraft 
performance  is  the  energy  state  model  developed  by  Bryson,  Desai,  and 
Hoffman  [37],  where  altitude  and  thrust  are  control  variables.  For 
many  aircraft,  this  energy  state  space  produces  a  hodograph  which  is 
not  convex.  The  physical  explanation  for  this  is  that  while  maintain¬ 
ing  a  constant  energy,  there  is  one  altitude  where  the  aircraft  flies 
aerodynamical ly  efficient  and  another  altitude  where  it  flies  thrust 
(or  power)  efficient.  This  dichotomy  leads  to  the  chattering  (or 
relaxed)  controls  solution  derived  by  Gilbert  [35] . 

To  eliminate  the  non-convexity  in  the  hodograph  and  preserve 
an  "optimal”  steady-state  cruise,  Schultz  and  Zagalsky  [6]  revised  the 
energy  state  model  so  that  altitude  becomes  a  state  variable  while 
flight  path  angle  and  thrust  form  the  control  space.  For  this  model 
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the  first  order  conditions  are  satisfied  using  controls  lying  in  a 
region  interior  to  the  control  space  where  the  legendre-Clebsch  con¬ 
dition  is  met  only  in  weak  form.  This  is  the  result  of  the 
Hamiltonian  being  a  linear  function  of  both  control  variables  and  is 
called  a  doubly,  singular  arc.  However,  Speyer  [7]  shows  that  the 
generalized  Legendre-Clebsch  (or  Kelley)  condition  [38]  is  not  satis¬ 
fied  for  this  case,  and  therefore,  the  steady-state  cruise  arc  is 
non-optimal . 

Schultz  [8]  once  more  revised  the  model;  this  time  so  the 
flight  path  angle  becomes  a  state  variable,  and  the  angle-of-attack 
and  thrust  form  the  control  space.  Now  with  this  model,  the  Kelley 
condition  is  satisfied.  Hence,  the  performance  index  with  respect 
to  the  controls  in  the  vicinity  of  the  steady-state  cruise  point 
appears  to  be  convex.  However,  once  again  Speyer  [5]  shows  that  the 
static  aircraft  cruise  is  non-minimizing.  Since  the  cruise  arc  in 
question  is  time  invariant,  a  transformation  to  the  frequency  domain 
can  be  made  which  greatly  simplifies  the  Jacobi  test,  the  final  and 
most  difficult  condition  to  apply.  Assuming  constant  mass,  infinite 
time,  and  fixed  end-points,  the  frequency  test  [5,  11,  18]  shows  that 
cyclic  solutions  produce  better  performance  than  the  steady-state 
solution  over  the  same  period.  However,  the  test  does  not  give  any 
indication  of  what  the  optimal  path  might  be  and,  as  yet,  this  remains 
an  interesting  and  unsolved  problem  [39]. 

Continuing  Speyer's  work,  Walker  [9]  searched  for  an  optimal 
periodic  flight  path.  However,  the  standard  numerical  optimization 
techniques  he  uses,  which  include  steepest  ascent  and  conjugate 
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gradient  methods,  do  not  converge  to  a  solution.  During  the  initial 
research  for  this  dissertation,  a  similar  effort,  using  a  hypersonic 
aircraft  dynamic  model,  also  failed  to  converge  to  a  solution.  This 
experience,  as  well  as  the  lack  of  suitable  illustrative  examples  in 
the  literature,  motivated  the  construction  of  a  less  complex  problem 
for  which  optimal  periodic  solutions  could  be  obtained  more  easily. 
The  assumption  made  here  is  that  the  characteristics  of  these  optimal 
periodic  processes  are  representative  of  those  for  the  general  class 
of  optimal  periodic  control  problems.  The  objective  ultimately  is  to 
develop  techniques  suitable  for  solving  problems  such  as  the  optimal 
cruise  trajectory  of  an  aircraft. 

5.2  Sample  Problem  Description 


A  particular  problem  has  been  constructed  to  represent  the 
general  class  of  optimal  periodic  control  problems  with  unrestricted 
period  as  defined  in  section  1.2.  The  sample  problem  is  described  in 
this  section  and  local  necessary  conditions  are  formulated  in  terms 
of  the  parameters  of  the  problem.  The  steady-state  solution  of  the 
problem  is  examined  for  optimality  using  the  frequency  test. 

A  statement  of  the  sample  problem  follows.  Find  the  period, 
T,  the  scalar  control,  u(*)>  and  the  state,  x(*)>  that  minimize  the 
performance  index, 


(5.1) 


subject  to  the  second  order  dynamic  constraints, 


■’-.i *Ub*tn> 
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^£1  (5.2) 

x2  =  U,  (5.3) 

and  to  the  periodicity  condition  applied  at  the  boundary, 

xx(T)  =  Xj(0) ,  (5.4) 

x2(T)  =  x2(0).  (5.5) 

The  performance  index  (5.1)  was  constructed  so  that  for 
particular  conditions,  periodic  solutions  would  exist  that  provided 
better  performance  than  the  steady-state  solution.  The  negative 
quadratic  term  provides  the  non-convexity  in  the  performance  criterion 
required  to  permit  solutions  other  than  the  static  solution.  The 
quartic  term  dominates  the  quadratic  term  for  large  excursions  of  x2, 
thereby  bounding  the  minimizing  solutions.  The  scalar  parameter,  b, 
weighting  the  control,  determines  the  nature  of  the  local  optimal 
solutions  as  will  be  shown  when  the  frequency  test  is  applied. 

First  the  variational  Hamiltonian  for  this  problem,  defined 
in  general  by  equation  (2.2),  may  be  written  as 
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*2u, 


(5.6) 


where  and  A2  are  the  Lagrange  multipliers  associated  with  the 
dynamic  constraints  (5.2)  and  (5.3).  Reducing  the  Euler-Lagrange 
equations,  (2.9)  through  (2.11),  by  solving  for  and  eliminating  the 


control,  u,  results  in  the  following  set  of  nonlinear  differential 


equations, 


x2  -  x2  -  Ar 


(5.7) 


(5.8) 


(5.9) 


(5.10) 


The  periodicity  conditions,  (5.4)  and  (5.5),  and  the  transversality 
condition  (2.13),  provide  the  boundary  conditions,  rewritten  below, 
for  this  two  point  boundary  value  problem. 


Xj(T) 


x2(T) 


Xj(T) 
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Xj  (0) , 


x2(0), 


V°), 


A2(0). 


(5.11) 


(5.12) 


(5.13) 


(5.14) 


The  form  of  the  optimal  periodic  control  is  given  by 


(5.15) 


where  A.”  satisfies  the  above  two-point  boundary  value  problem.  All 
solutions  to  this  problem  are  extrema  of  the  sample  optimal  periodic 
control  problem.  Also  locally  optimizing  solutions  of  the  problem 
must  satisfy  the  additional  transversality  condition  (2.14), 


3L-*  *yW' 
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This  condition  extracts  from  the  previously  obtained  set  of  all 
extrema  a  subset  of  which  each  element  has  been  extremized  with 
respect  to  an  unconstrained  period. 

Now  consider  the  second  order  necessary  conditions.  By 
restricting  the  control  weighting  parameter  such  that  b  >  0,  the 
Legendre-Clebsch  condition  is  always  satisfied  in  its  strong  form, 
i.e.,  Huu  =  b  >  0.  The  Weierstrass  condition  is  also  always  satisfied 
since  the  Hamiltonian,  (5.6),  is  regular,  i.e.,  it  has  only  one  minimum 
with  respect  to  u. 

For  the  sample  problem,  there  is  only  one  steady  state  solu¬ 
tion  which  satisfies  all  of  the  above  necessary  conditions  for  an 
optimal  solution, 

Xj  =  x2  =  u  =  0.  (5.17) 

By  using  the  frequency  test  [11],  the  range  of  the  parameter,  b,  for 
which  this  static  path  is  actually  minimizing  can  be  determined.  The 
matrix  expressions,  corresponding  to  those  in  (2.22)  of  the  frequency 
test,  in  terms  of  the  sample  problem  are 


(5.18) 


where  f  is  the  functional  representation  of  the  dynamic  constraints, 
(5.2)  and  (5.3).  The  Hermitian  matrix  (2.23),  a  form  of  the  second 
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variation  transformed  into  the  frequency  domain,  becomes 


n(u») 


(5.19) 


after  eliminating  terms  multiplied  by  the  zero  vector,  Hxu.  The 
matrix  expression  may  be  further  reduced  by  introducing  the  remaining 
relationships  from  (5.18)  into  (5.19),  which  gives 


n(u>)  =  ±  -  i  +  b. 

U)  U) 


(5.20) 


As  indicated  in  section  2.4,  the  steady  state  solution  is 
minimizing  if  the  matrix  n(u>)  is  non-negative  definite  for  all  values 
of  the  frequency,  u>.  The  minimum  value  for  n(u>)  in  (5.20)  occurs  when 
u»  =  J2.  As  a  result,  when  b  Z  fc,  the  matrix  n(ui)  is  non-negative 
definite,  and  the  static  solution  is  minimizing  for  this  range  of  the 
parameter,  b.  However,  in  the  range  0  £  b  <  there  are  frequencies 
for  which  the  matrix  is  negative.  For  these  values  of  b,  there  are 
periodic  solutions  which  improve  performance.  The  extraction  and  char¬ 
acterization  of  these  solutions  are  important  objectives  of  this  work. 
Therefore,  a  thorough,  numerical  and  analytical  study  of  the  sample 
problem,  formulated  in  this  section,  is  conducted  in  the  remaining 
chapters  of  this  work.  But  first,  some  special  relationships  result¬ 
ing  from  the  symmetry  and  order  of  the  sample  problem  are  developed  in 
the  last  two  sections  of  this  chapter. 

5.3  The  Trace  of  the  Monodromy  Matrix 

The  trace  of  a  matrix  is  an  easily  obtained  and  a  quite 


powerful  indicator  of  some  intrinsic  properties  of  a  matrix.  Because 
of  the  relationship  of  the  trace  and  the  eigenvalues  of  a  matrix,  the 
stability  of  a  system  associated  with  the  eigenvalues  of  the  matrix 
characterizing  the  system  is  frequently  determined  by  obtaining  its 
trace.  As  was  derived  in  section  4.2  for  the  general  optimal  periodic 
control  problem,  locally  optimizing  solutions,  except  for  isolated 
cases,  belong  to  the  set  of  unstable,  periodic  solutions  to  the  Euler- 
Lagrange  equations,  (2.9)  through  (2.11).  For  the  case  of  a  fourth 
order,  symmetric  system,  this  result  is  demonstrated  by  a  very  simple 
relationship  involving  the  trace  of  its  monodromy  matrix  and  elements 
of  the  Riccati  variable  associated  with  the  system. 

The  sample  problem  formulated  in  this  chapter  generates  a 
fourth  order  system  of  differential  equations,  (5.7)  through  (5.10). 
The  system  is  symmetric,  satisfying  the  conditions  for  symmetry  given 
in  appendix  B  by  equation  (B.2).  The  monodromy  matrix  for  the  sample 
problem,  as  well  as  for  any  other  fourth  order,  symmetric  system,  may 
be  expressed  by, 

r  r 

r  n  12 

r  r 

21  22 

This  special  form  of  the  monodromy  matrix  for  a  symmetric  system  is 
derived  in  appendix  B. 

It  is  apparent  from  examining  the  monodromy  matrix  (5.21) 
that  only  its  first  two  columns  need  be  determined  to  calculate  the 
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trace  of  the  matrix, 

trace  P  =  2(a  +  d) .  (5.22) 

Using  the  relationship,  (3.48),  between  the  trace  and  the  eigenvalues 
of  a  matrix,  the  reciprocity  property  of  the  eigenvalues  of  the 
transition  matrix  (section  3.3),  and  the  property  that  at  least  two 
eigenvalues  of  the  monodromy  matrix  are  unity  (section  3.4)  permits 
expressing  the  trace  of  the  monodromy  matrix  (5.22)  in  terms  of  one 
unknown  eigenvalue,  p, 

2(a  +  d)  =  1  +  1  +  p  +  .  (5.23) 

From  the  similarity  transformation  in  Theorem  4.1,  another 
expression  for  the  trace  of  the  monodromy  matrix  may  be  written  as 
follows, 

trace  T  =  trace  +  P^pJ  +  trace  |V22  -  Prj2J,  (5.24) 

where  P  is  the  initial  value  matrix  of  the  Riccati  variable.  It  was 
shown  in  section  4.2  that  this  similarity  transformation  also  parti¬ 
tions  the  eigenvalues  of  P  into  the  two  submatrices  on  the  diagonal 
of  the  transformed  matrix.  The  eigenvalues  of  the  submatrix, 

[r^  +  ri2P]  ,  are  the  reciprocal  of  the  eigenvalues  of  the  other  sub¬ 
matrix,  [r22  -  Prj2l.  As  a  result,  an  expression  for  the  trace  of  the 
submatrix  (P^  +  ri2P]  may  written  as 

trace  jYn  +  P^pj  =  1  +  M-  (5.25) 
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Since  the  Riccati  matrix  is  symmetric,  it  can  be  expressed 
in  terms  of  its  scalar  elements  as  follows, 


P 


(5.26) 


Using  equations  (5.24)  and  (5.26)  to  perform  the  matrix  operations 
indicated  by  equation  (5.25)  gives  the  result, 


a  +  d  +  eP1  +  gP3  =  1  +  p.  (5.27) 

Subtracting  equation  (5.27)  from  (5.23)  gives 

a  +  d  -  ePj  -  gP3  =  1  +  J  .  (5.28) 

Eliminating  the  eigenvalue  \i  from  the  two  equations  (5.27)  and  (5.28) 
and  then  rearranging  the  resulting  expression  provides  the  following 
useful  relationship, 

(a  +  d  -  l)2  -  (ePj  +  gP3)2  =  1.  (5.29) 

This  equation  clearly  isolates  the  range  of  values  that  the  trace  of 
the  monodromy  matrix  may  assume  when  the  elements  of  the  Riccati  matrix 
are  real.  This  leads  to  a  necessary  condition  for  local  optimality 
of  extremal  solutions  to  the  optimal  periodic  control  problem  which 
may  be  expressed  in  terms  of  elements  of  the  trace  of  the  monodromy 
matrix, 


a  +  d  £  2 . 


(5.30) 


The  necessary  condition  (5.30)  for  fourth  order,  symmetric  systems 
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excludes  any  solution  whose  eigenvalues  lie  on  the  unit  circle  at 
points  other  than  the  critical  points,  (1,1)  or  (-1,-1).  This  is  in 
full  agreement  with  the  more  general  results  derived  in  section  4.2. 

5 . 4  A  Solution  to  the  Riccati  Equation 

Another  important  condition  for  optimality  involves  the  ex- 
istance  of  the  Riccati  variable  over  the  full  time  duration  spanned  by 
the  optimal  solution.  When  the  solution  is  periodic,  only  an  inter¬ 
val  of  one  period  need  be  considered.  In  this  section  the  algebraic 
Riccati  equation,  derived  in  section  4.1,  is  evaluated  for  fourth 
order,  symmetric  systems.  Two  initial  value  matrices  for  the  Riccati 
variable  are  derived  in  terms  of  elements  of  the  monodromy  matrix. 
These  results  explicitly  determine  the  only  periodic  solutions  to  the 
Riccati  matrix  differential  equation  (2.15).  With  the  initial  condi¬ 
tions,  the  existance  of  both  sets  of  Riccati  variables  over  a  full 
period  may  be  examined  by  integrating  the  differential  equation. 

The  algebraic  Riccati  equation  (4.16)  provides  the  general 
relationship  between  the  elements  of  the  monodromy  matrix  and  the 
initial  conditions  for  periodic  solutions  to  the  Riccati  differential 
equation.  This  relationship  may  be  specialized  for  fourth  order, 
symmetric  systems  by  expanding  the  algebraic  equation,  using  the  par¬ 
ticular  expressions  for  the  monodromy  matrix  and  the  Riccati  matrix 
given  by  equations  (5.21)  and  (5.26).  The  results  of  this  matrix 
equation  expansion  are  expressed  below  in  the  equivalent  form  of  four 
simultaneous,  algebraic  scalar  equations, 
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eP2  +  gP2  +  2cP2  -  h  =  0,  (5.31) 

-£P2  *  eP>P2  *  fPlP3  *  *P2P3 

+  (d  -  a)P2  +  bP1  +  cP3  -  i  =  0,  (5.32) 

2 

fp  +  eP  P  -  fP  P  +  eP  P 
2  1  2  1  3  8  2  3 

+  (a  -  d)P2  +  bPj  +  cP3  +  i  =  0,  (5.33) 

eP2  +  gP2  +  2bP2  -  j  =  0.  (5.34) 

Because  the  Riccati  matrix  is  symmetric,  there  are  three  unknown  ele¬ 
ments  in  the  matrix.  Therefore,  only  three  of  the  above  equations  are 
independent. 

A  particular  combination  of  these  equations  results  in  an 
explicit  solution  for  the  initial  condition  of  the  Riccati  variable, 
P2,  expressed  in  terms  of  elements  of  the  monodromy  matrix  only. 
Subtracting  equation  (5.33)  from  equation  (5.32)  and  multiplying  the 
result  by  the  monodromy  matrix  elements  e  and  g  produces  the  follow¬ 
ing  equation, 

f(2egP1P3)  -  2efgP2  +  2eg(d  -  a)P2  -  2egi  =  0.  (5.35) 

Multiplying  equation  (5.31)  by  the  elements  e  and  f  and  adding  the 
result  to  that  of  multiplying  equation  (5.34)  by  the  elements  f  and 
g  gives  a  second  equation, 

f(e2P^  +  g2P2  )  +  2efgP2  +  2f(ce  +  bg)P2 

-  f (eh  +  gj)  =  0.  (5.36) 
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The  unknowns,  P^  and  P^,  can  be  eliminated  by  adding  the  two  equations, 
(S.3S)  and  (5.36),  and  replacing  the  term  involving  P^  and  P^  in  the 
result  with  its  equivalent  expression, 

(eP1  +  gP3)2  =  (a  +  d  -  l)2  -  1,  (5.37) 


derived  in  the  previous  section.  This  results  in  a  simple,  algebraic 
expression  for  P2  involving  only  the  elements  of  the  monodromy  matrix, 


p  _  2egi  +  f[eh  +  gj  +  1  -  (a+d-1)2] 
2  “  2  [f(ce  +  bg)  +  eg(d  -  a)] 


(5.38) 


It  is  interesting  that  the  initial  condition  for  this  component  of 
the  Riccati  matrix  is  single-valued.  This  implies  that  there  is  only 
one  starting  condition,  P^,  associated  with  .each  initial  condition, 
tQ,  for  the  optimal  periodic  control  problem,  that  results  in  a 
periodic  solution  to  the  corresponding  Riccati  differential  equation. 

Expressions  for  the  initial  conditions  P^  and  P^  can  now  be 
obtained  by  using  the  value  for  P^  from  equation  (5.38)  in  equations 
(5.31)  and  (5.34),  rearranged  below, 


P 


1 


P 


3 


Jh  -  2cP2-  gP*j*  _  ,nd 

■ 


(5.39) 


(5.40) 


The  indeterminacy  of  the  corresponding  signs  in  equations  (5.39)  and 
(5.40)  may  be  alleviated  by  the  following  expression  relating  Pj  and 
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eP  +  b 

p3  =  --T—  pr  <5-41> 

gP2  +  c 

This  relationship  was  obtained  by  adding  equations  (5.32)  and  (5.33) 
and  rearranging  the  result. 

The  expressions  (5.38)  through  (5.41)  provide  all  of  the 
initial  conditions  which  result  in  a  periodic  solution  to  the  Riccati 
differential  equation,  (2.15).  In  general,  and  at  most,  there  are 
two  sets  of  these  starting  conditions.  Corresponding  to  each  set  is 
a  single  set  of  periodic  solutions  which  satisfy  the  differential 
equation.  It  may  be  deduced  from  Theorem  4.1  that  the  two  periodic 
matrix  solutions  have  properties  that  are  related  to  the  pairing  of 
the  eigenvalues  in  the  submatrices  on  the  diagonal  of  the  transformed 
matrix  (4.19).  These  results  are  explicitly  shown  in  the  numerical 
investigation  of  the  next  chapter. 


CHAPTER  6 


NUMERICAL  INVESTIGATION 

A  numerical  investigation  of  the  sample  optimal  periodic 
control  problem  illustrates  the  extraordinary  complexity  of  this  class 
of  optimization  problem.  Periodic  solutions  to  the  Euler-Lagrange 
equations  are  obtained  in  the  first  section  that  are  continuously 
related  forming  a  one-parameter  family  of  solutions.  In  section  two, 
particular  solutions  are  identified  as  bifurcation  points  by  a  measure 
of  their  stability.  Additional  families  of  solutions  are  obtained  that 
branch  from  the  original  family  at  the  bifurcation  points,  and  charac¬ 
teristics  that  distinguish  their  solutions  are  examined.  Local  minimal 
solutions  are  determined  in  the  next  section.  The  performance  index 
associated  with  each  periodic  extremum  is  computed  and  the  necessary 
condition  for  optimization  with  respect  to  the  period  is  applied.  In 
section  four,  the  sufficiency  condition  associated  with  the  existence 
of  the  Riccati  variable  is  examined.  Finally,  in  the  last  section 
the  periodic  regulator  which  was  developed  in  chapter  four  is  demon¬ 
strated  using  the  local  optimal  periodic  solution  associated  with  the 
principal  family. 

6. 1  Principal  Family  of  Periodic  Solutions 

The  Euler-Lagrange  equations  in  the  form  of  a  two-point 
boundary  value  problem,  equations  (5.7)  through  (5.14),  are  solved  by 
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numerical  techniques  in  this  section.  This  problem  was  constructed 
to  be  as  simple  as  possible  and  still  produce  periodic  solutions. 
However,  the  resulting  solution  to  the  problem  is  remarkably  complex. 
An  infinity  of  periodic  solutions  is  shown  to  exist  for  the  problem 
forming  one-parameter  families.  A  distinguishing  characteristic  of 
solutions  in  a  family  is  the  number  of  positive  crossings  of  the  axis 
in  one  period.  The  principal  family,  which  originates  at  the  static 
equilibrium  solution  of  the  problem,  is  determined  and  some  character¬ 
istics  of  its  solutions  are  examined. 

An  existing  computer  program,  utilizing  a  "shooting  method", 
was  used  to  compute  the  periodic  solutions  obtained  in  this  work. 

The  program  originally  was  developed  by  Roger  A.  Broucke  in  his  study 
of  periodic  orbits  near  equilibrium  points  in  the  restricted  problem 
of  three  bodies  [24].  It  was  designed  to  find  periodic  orbits 
(solutions)  for  any  fourth  order  Hamiltonian  system,  given  a  set  of 
initial  conditions  within  the  range  of  convergence  for  the  system. 
Suitable  initial  conditions  may  be  found  using  one  of  several  search 
routines  available  in  the  program.  After  determining  a  periodic 
solution,  the  program  uses  an  interpolation  routine  to  predict  initial 
conditions  for  a  neighboring  periodic  solution.  In  this  manner  it  de¬ 
termines  and  follows  a  family  of  periodic  solutions.  The  scheme  is 
analogous  to  the  principal  of  analytic  continuation,  described  by 
Szebehely  [23],  and  it  exploits  the  inherent  properties  of  the 
Hamiltonian  system  and  its  associated  monodromy  matrix. 

The  periodic  solutions  obtained  for  this  problem  are  unique¬ 
ly  identified  by  their  initial  conditions  as  was  expressed  previously 
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in  chapter  three.  This  provides  a  very  useful  means  of  categorizing 
results  and,  in  the  case  of  fourth  order  systems,  illustrating  them. 

In  figures  6.1a  and  6.1b,  plots  of  the  principal  family  of  periodic 
solutions  are  shown.  Each  point  of  the  plots  represents  the  initial 
conditions  of  a  periodic  solution  to  the  Euler-Lagrange  equations  of 
the  sample  problem.  In  this  case,  the  initial  conditions  for  the 
variables  and  Aj,  are  arbitrarily  chosen  to  be  zero.  This  results 
from  the  order  of  the  system  of  differential  equations  being  reducible 
by  two  as  expressed  in  chapter  three.  The  remaining  two  initial 
conditions  are  x^,  the  abscissa  of  both  graphs,  and  implied  by  the 
Hamiltonian  in  figure  6.1a  and  explicitly  identified  on  the  ordinate 
in  figure  6.1b.  A  control  weighting  parameter  of  b  =  0.1  has  been 
used  in  all  examples  of  this  chapter. 

The  period  of  the  solutions  comprising  the  principal  family 
ranges  continuously  from  approximately  2.11  to  5.92.  The  two  extrema 
correspond  to  the  fast  and  slow  frequencies  represented  by  the  eigen¬ 
values  of  the  Euler-Lagrange  equations  linearized  about  the  static 
equilibrium  solution.  This  relationship  to  the  period  is  shown  for  x^ 
in  figure  6.2a  and  for  H  in  figure  6.2b.  Traversing  the  right  half 
plot  of  figure  6.1a  in  a  counter-clockwise  direction,  beginning  at 
the  origin,  corresponds  to  solution  with  increasing  period  in  figures 
6.2a  and  6.2b. 

6.2  Bifurcation  Points  and  New  Families 

The  stability  of  a  periodic  solution  is  an  important  inher¬ 
ent  characteristic  of  the  solution.  It  is  also  a  key  element  in  the 
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relationship  between  families  of  solutions.  In  this  chapter,  as  in 
previous  chapters,  a  periodic  solution  is  stable  if  no  eigenvalue  of 
its  monodromy  matrix  lies  outside  the  unit  circle.  Otherwise,  it  is 
unstable.  A  measure  of  the  stability  of  a  periodic  solution  is  the 
stability  index,  k,  which  is  defined  as  the  trace  of  the  monodromy 
matrix  of  the  solution  less  two,  for  a  fourth  order  system.  An  equiv¬ 
alent  definition  is  k  equals  the  sum  of  the  remaining  eigenvalues  of 
the  monodromy  matrix  after  eliminating  the  two  unity  eigenvalues. 

The  relationship  of  the  stability  index  to  the  stability 
of  a  periodic  solution  for  a  fourth  order  system  is  identified  by  the 
following  inequalities: 


k  g  2, 

stable  solution; 

(6.1) 

k  >  2, 

unstable  solution. 

(6.2) 

It  is  interesting  to  note  that  the  stability  index,  k,  is  precisely 
that  quantity,  a  +  d,  which  occurs  in  the  inequality  equation  (5.30) 
expressing  a  necessary  condition  for  optimality  of  a  symmetric,  fourth 
order  periodic  control  system. 

An  important  relationship  used  by  Henon  [27]  and  Contopoulos 
[26]  expresses  properties  of  the  solutions  of  a  family  in  terms  of  a 
critically  stable  solution,  i.e.,  all  eigenvalues  equal  one.  The  re¬ 
lationship  may  be  written  in  terms  of  critical  values  of  the  stability 
index, 

kc  =  2  cos(2"n),  (6.3) 
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where  m  and  n  are  both  integers  associated  with  properties  of  periodic 
solutions  with  stability  index  equal  to  kc>  One  importance  of  this 
relationship  is  that  solutions  of  a  family  which  may  be  common  to 
other  branching  families  can  be  identified  by  their  stability  index. 
These  common  solutions  are  frequently  called  bifurcation  points. 

Some  points  for  the  principal  family  are  plotted  in  figure  6.3  and 
are  identified  by  the  notation  n/m  Fi,  where  n  and  m  are  the  integers 
above  (/m  is  omitted  when  m  is  one),  F  is  the  family  designator,  and 
i  distinguishes  equal  critical  values.  The  stable  and  unstable 
solutions  of  the  family  are  identified  by  a  broken  line  and  a  solid 
line  respectively.  The  unstable  solutions  have  no  bifurcation  points, 
according  to  equation  (6.3),  and  they  represent  extrema  of  the  problem 
which  also  satisfy  the  necessary  condition  for  optimality,  equation 
(5.30). 

Branching  families  are  found  by  searching  for  new  periodic 
solutions  in  the  vicinity  of  bifurcation  points.  The  new  families  are 
then  traced  in  the  same  manner  as  the  principal  family.  A  few  new 
families  branching  from  the  principal  family  are  shown  in  figure  6.4. 
Stable  and  unstable  solutions  are  again  identified  by  broken  and  solid 
lines  respectively.  Only  the  lower  half  of  the  plots  are  shown  because 
this  region  provides  the  solutions  of  most  interest  to  the  optimal 
periodic  control  problem.  Solutions  that  satisfy  the  necessary 
condition  for  optimal  period,  H  =  J°,  and  improve  performance  com¬ 
pared  to  the  steady  state  must  lie  in  the  region  H  <  0.  New  families 
are  identified  using  their  generating  bifurcation  point  identifier 
preceded  by  an  F. 
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Two  additional  levels  of  branching  are  shown  on  the  plot  in 
figure  6.5,  which  is  a  detailed  enlargement  of  a  section  of  the 
previous  plot.  This  process  may  be  continued  indefinitely,  limited 
only  by  the  finite  word  length  of  the  computer  and  the  persistence 
of  the  investigator.  It  is  quite  apparent  that  large  regions  of  the 
phase  space  are  densely  packed  with  families  of  periodic  solutions, 
each  an  extremum  of  the  optimal  periodic  control  problem. 

Before  proceeding  to  the  next  section  and  further  investigat¬ 
ing  the  sample  problem,  some  important  characteristics  of  the  periodic 
solutions  of  a  family  are  examined.  A  characteristic  that  distin¬ 
guishes  one  family  from  another  is  the  number  of  axis  crossings  that 
occur  in  the  same  direction  during  one  period  of  a  solution.  This  is 
frequently  expressed  as  the  number  of  arcs  of  a  solution.  Referring 
back  to  the  equation  (6.3),  the  integer,  n,  is  the  ratio  of  the  arcs 
of  a  solution  on  the  branching  family  to  the  arcs  of  a  solution  on  the 
principal  (or  branched  from)  family. 

This  characteristic  is  easily  discernible  in  the  following 
sequence  of  time  history  plots  and  phase  plots  of  solutions  represent¬ 
ing  different  families.  In  figure  6.6,  a  time  history  plot  of  the 
variable  x^  for  three  solutions  on  the  principal  family  are  shown. 

The  plots  a,  b,  and  c  correspond  to  solutions  at  the  bifurcation 
points  2A,  5/2A,  and  1A.  This  sequence  illustrates  the  variation  of 
Xj  in  amplitude  and  period  along  the  principal  family  FA.  In  figure 
6.7,  a  similar  sequence  of  plots  is  shown  for  the  branching  family 
F2A.  The  evolution  of  a  two  arc  solution  from  two  periods  of  a  one 
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arc  solution  beginning  at  the  bifurcation  point  2A  is  identifiable  in 
the  sequence  of  plots  a,  b,  and  c.  The  same  sequence  of  plots  for  the 
family  FS/2A  is  depicted  in  figure  6.8.  Plot  a,  corresponds  to  five 
periods  of  the  one  arc  solution  at  the  bifurcation  point  5/2A.  Plot 
c,  corresponds  to  the  solution  at  the  bifurcation  point  2(5/2A)  as 
shown  on  the  branching  family  identified  in  the  detailed  enlargement, 
figure  6.5.  Plot  b,  is  an  intermediate  solution  on  the  family.  The 
final  time  history  sequence  of  plots,  figure  6.9,  shows  the  family 
F2(5/2A).  Two  periods  of  the  five  arc  solution  at  the  bifurcation 
point  2(5/2A)  are  shown  in  plot  a.  The  remaining  two  plots,  b  and  c, 
show  the  evolution  of  a  ten  arc  solution. 

A  useful  and  probably  more  frequently  used  plot  of  periodic 
solutions  is  one  providing  phase  space  relationships.  In  figure  6.10, 
phase  plots  relating  the  variable  x^  to  x^,  corresponding  to  the 
families  FA,  F2A,  F5/2A,  and  F2(5/2A)  are  shown  in  plots  a  through  d 
respectively.  Each  of  the  four  plots  represent  the  corresponding 
solutions  in  plot  c,  of  the  previous  time  history  plots. 

Considerable  additional  work  has  been  completed  by  Henon  [27] 
and  Contopoulos  [26]  classifying  bifurcation  points.  For  certain 
conditions,  which  exist  in  this  problem,  trifurcations  are  predicted. 
Frequently,  one  of  the  two  or  three  new  families  is  composed  of  non- 
symmetric  solutions.  For  such  a  case  the  family  does  not  lie  in  the 
symmetric  plane  (x^,  0,  0,  •  Further  classification  of  critical 

solutions  of  the  first  kind  (kc  =  -*-2)  and  the  second  kind  (k£  =  -2) 
are  summarized  in  appendix  C. 
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6.3  Optimized  Solutions  for  Free  Period 

The  necessary  condition  for  an  optimal  periodic  solution  when 
its  period  is  unconstrained  is  applied  in  this  section  to  the  extremal 
solutions  obtained  in  the  previous  sections.  This  condition  screens 
the  infinity  of  solutions  forming  families  so  that  only  individual 
solutions  of  the  families  remain.  Some  interesting  observations  from 
the  results  obtained  that  pertain  to  the  performance  index  are  expres¬ 
sed  concluding  the  section. 

The  performance  index  is  calculated  for  solutions  of  the 
principal  family  and  the  results  are  plotted  in  figure  6.11.  For 
these  plots  the  solid  line  represents  the  initial  conditions  of 
periodic  solutions  to  the  Euler-Lagrange  equations  as  in  figure  6.1a, 
where  the  initial  conditions,  x2  and  Xj,  are  zero,  the  initial  con¬ 
dition  Xj  is  directly  determined  on  the  abscissa,  and  the  initial 
condition  is  implied  by  H  on  the  ordinate.  The  broken  line  identi¬ 
fies  on  the  ordinate  the  performance  index,  associated  with  each 
periodic  solution  corresponding  to  the  initial  conditions  represented 
by  Xj.  At  the  intersection  of  the  two  plots  is  the  solution  for  which 
the  condition  H  =  J°  is  satisfied.  The  upper  half  of  the  plots,  H  >  0, 
has  been  omitted  to  simplify  the  figure.  Similarly,  results  for  faiaily 
F1A  are  plotted  in  figure  6.12.  The  expanded  ordinate  in  figure  6.12b 
better  illustrates  the  extrema.  The  two  apparent  local  minima  (ex¬ 
cluding  those  resulting  from  the  symmetry  of  the  problem)  for  this 
family  are  actually  the  same  solution  intersecting  the  surface  of 
initial  conditions  of  the  plot  in  two  points. 
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The  performance  index  as  a  function  of  the  period  is  shown 
in  figure  6.13  on  plot  a  for  the  principal  family  FA  and  on  plot  b  for 
the  branching  family  F2A.  The  extrema  of  the  performance  index ,  deter¬ 
mined  from  the  plots  in  figure  6.13,  agree  with  the  corresponding 
results  in  figures  6.11  and  6.12. 

Periodic  solutions  satisfying  all  the  first  order  necessary 
conditions  for  a  local  optimum  with  unconstrained  period  are  plotted 
for  the  principal  family  FA  and  the  branch  families  F2A,  F3A,  and 
F5/2A.  In  figure  6.14  phase  plots  relating  the  variable  x^  to  x^  are 
shown  for  the  four  local  minimum  solutions.  A  time  history  of  the 
control  variable  corresponding  to  each  of  the  previous  four  local  min¬ 
imum  solutions  is  plotted  in  figure  6.15.  The  initial  conditions, 
optimal  period,  and  performance  index  for  the  local  minima  of  the 
four  families  are  summarized  in  table  6.1. 

Concluding  this  section  are  three  interesting  observations 
made  during  the  numerical  investigation  which  may  or  may  not  be  true 
in  general.  First,  every  periodic  solution  to  the  Euler-Lagrange 
equations  obtained  in  this  investigation  provided  an  improvement  when 
compared  to  the  performance  index  of  the  steady-state  solution. 

Second,  the  performance  index  is  a  relative  maximum  for  the  branch 
family  when  the  Hamiltonian  is  also  an  extremum  with  respect  to  the 
initial  condition.  Finally,  there  are  no  solutions  with  period  less 
than  that  associated  with  the  fast  frequency  eigenvalue  (from  linear 
analysis  about  the  static  equilibrium  point)  that  improves  the  steady- 
state  performance. 
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6.4  Sufficiency  Condition 

The  solutions  of  the  previous  section  must  complete  one  final 
check  for  optimality.  The  Riccati  variable  associated  with  each  solu¬ 
tion  must  exist  over  a  complete  period.  The  existence  of  the  Riccati 
variable  is  examined  by  two  equivalent  methods  which  cong>letes  the 
test  of  the  sufficiency  conditions  for  optimality  of  the  periodic 
control  problem  expressed  by  Theorem  4.2. 

The  initial  conditions  for  each  periodic  solution  of  the 
Riccati  differential  equation  (2.15)  are  determined  directly  from  the 
algebraic  expressions  (5.38)  through  (5.41).  Two  solutions  exist  for 
fourth  order  symmetric  systems.  They  are  obtained  by  integrating  the 
Riccati  equation  over  one  period  using  the  two  sets  of  starting  condi¬ 
tions  just  obtained.  An  equivalent  procedure  involves  evaluating  the 
expression  f4>11(t,0)  +  4>12(t,0)P(0)]  from  equation  (4.13),  and  test¬ 
ing  for  a  zero  value.  The  existence  of  the  Riccati  variable  is 
implied  by  no  zero  values  in  a  one  period  interval. 

The  Riccati  variable  corresponding  to  the  optimum  periodic 
solution  of  the  principal  family  was  found  to  exist  over  the  one  period 
interval,  thus  confirming  the  local  optimality  of  the  solution.  The 
time  histories  of  the  three  elements  comprising  the  Ricatti  variable 
are  plotted  for  both  sets  of  initial  conditions  in  figures  6.16  and 
6.17. 

Each  of  the  remaining  extremal  solutions  from  the  previous 
section  evolved  from  a  branch  family.  At  some  point  in  the  one  period 
interval  of  each  of  these  extremal  solutions,  the  corresponding  Riccati 
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variable  escaped  (did  not  exist).  Therefore,  these  solutions  are  not 
locally  optimizing.  As  a  result,  out  of  the  many  extremal  solutions 
that  have  been  found  for  this  problem,  only  the  one  solution  on  the 
principal  family  satisfies  each  condition  for  local  optimality  of  a 
periodic  solution  with  unspecified  period. 

The  existence  of  the  Riccati  variable  was  also  examined  for 
other  extrema  in  the  unstable  region  (stability  index,  k  >  2)  of  the 
principal  family  and  the  branch  families.  The  basic  results  were  the 
same.  Corresponding  to  extrema  of  the  principal  family,  the  Riccati 
variable  was  found  to  exist  during  the  entire  one  period  interval  for 
the  sampling  of  solutions  tested.  However,  corresponding  to  extrema 
of  the  branch  families,  the  Riccati  variable  escaped  at  some  time  in 
the  interval  for  each  of  the  solutions  tested. 

These  results  indicate  that  the  extremal  solutions  in  the  un¬ 
stable  region  of  the  principal  family  are  local  optimal  solutions  with 
respect  to  a  specified  period.  The  range  of  periods  corresponding  to 
the  optimal  solutions  extend  from  approximately  2.67  to  4.99.  Optimal 
solutions  for  specified  periods  outside  this  range  are  as  yet  unknown. 
If  it  could  be  assumed  that  n  cycles  of  an  optimal  periodic  solution 
with  respect  to  a  specified  period,  T,  were  also  optimal  for  a  speci¬ 
fied  period,  nT,  there  would  still  be  gaps  in  the  range  of  periods  for 
which  optimal  solutions  exist.  It  is  strongly  suspected  that  nonsym- 
metric  solutions,  those  whose  initial  conditions  lie  outside  the  sym¬ 
metric  surface  investigated  in  this  work,  will  provide  additional 
local  optimal  solutions  for  the  problem. 
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6.5  Periodic  Regulator 

A  periodic  regulator  was  developed  in  chapter  four  using  the 
neighboring  optimal  feedback  control  law,  (4.57).  In  this  section, 
the  regulator  is  implemented  using  the  optimal  periodic  solution  of 
the  previous  section  as  the  nominal  solution.  The  regulator  is  tested 
by  perturbing  an  initial  state  and  observing  the  results. 

The  control  diagram  for  the  periodic  regulator  is  shown  in 
figure  6.18a.  The  plant  (or  system)  for  this  problem  is  defined  by 
the  dynamic  equations  (5.2)  and  (5.3)  rewritten  below, 

Xj  =  x2,  (5.2) 

x2  =  u.  (5.3) 

The  gains  for  the  regulator  are  determined  from  the  optimal  neighbor¬ 
ing  feedback  control  law  (4.57)  rewritten  below, 

6u°  =  (H°  )_1(f°P  +  H°  )6x.  (4.57) 

uu  u  ux 

For  this  problem  the  gains  reduce  to 


G 


b^2  **3^  ’ 


(6.4) 


where  b  =  0.1  and  the  elements  of  the  Riccati  variable  are  chosen  from 
one  of  the  two  possible  sets  shown  in  figures  6.16  and  6.17.  One  of 
the  sets  causes  a  perturbed  system  to  converge  back  to  the  nominal 
path  in  a  limit  cycle;  the  other  causes  the  system  to  diverge. 

The  storage  block  of  the  regulator  is  composed  of  a  look-up 
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table  and  a  least  squares  error  measurement  algorithm.  The  optimal 
periodic  solution  from  the  previous  section  is  used  for  determining 
the  nominal  values  in  the  table.  The  terms  in  the  table  represent 
the  nominal  path  plotted  in  figure  6.14a,  the  nominal  control  plotted 
in  figure  6.15a,  and  the  gains  determined  above.  The  table  spans  a 
full  period  and  each  increment  represents  a  time  interval.  As  shown 
in  figure  6.18b,  the  error  measurement  algorithm  used  is  time  indepen¬ 
dent.  Given  the  actual  state  coordinates,  the  algorithm  determines 
from  the  table  the  nominal  path  coordinates  which  represent  the  short¬ 
est  distance  between  the  actual  state  and  the  nominal  path.  The 
nominal  control  and  the  gains  which  correspond  in  the  table  to  the 
nominal  path  coordinates  just  determined  are  used  in  this  increment 
of  the  control  cycle. 

The  regulator  is  demonstrated  by  initially  perturbing  the 
state  and  then  observing  the  output  of  the  plant.  Two  cases  are  pre¬ 
sented  for  comparative  purposes.  For  both  cases  the  output  of  the 
plant  (actual  path)  is  represented  by  a  broken  line  and  the  nominal 
path  a  solid  line. 

In  the  first  case  the  periodic  regulator  is  demonstrated  by 
closing  the  feedback  loop  and  using  the  set  of  divergent  gains. 
Cooresponding  to  these  gains  are  the  set  of  Riccati  variable  elements 
shown  in  figure  6.17.  The  initial  condition,  Xj(0),  is  perturbed  by 
0.01  per  cent.  The  result,  shown  in  figure  6.19a,  illustrates  the  di¬ 
vergent  character  of  this  set  of  gains. 

The  final  case  implements  the  periodic  regulator  with  its 
set  of  convergent  gains.  The  Riccati  variables,  shown  in  figure  6.16, 
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are  used  to  determine  these  gains.  For  this  case  the  initial  condition 
is  perturbed  from  the  nominal  by  50  percent.  The  strongly  convergent 
behavior  of  the  periodic  regulator  is  shown  in  figure  6.19b. 


d  . 


CHAPTER  7 


AN  ASYMPTOTIC  SOLUTION 

An  approximate  analytical  solution  to  the  example  problem, 
defined  in  chapter  five,  is  determined  in  this  chapter.  Using  a  per¬ 
turbation  technique,  the  solution  is  generated  in  the  form  of  an  asymp¬ 
totic  series  expansion.  Both  the  dependent  and  independent  variables 
are  expanded  about  a  small  perturbation  parameter.  The  success  of  this 
technique  is  very  much  dependent  on  the  selection  of  the  functional 
form  of  the  parameter.  In  the  first  section,  the  example  problem  is 
reformulated  in  terms  of  expanded  variables.  An  expression  in  the 
form  of  a  Fourier  series  expansion  for  the  solution  is  determined  in 
the  second  section.  Then  in  the  following  section,  the  conditions  for 
an  optimal  period  are  applied.  This  results  in  asymptotic  series 
expressions  for  the  optimal  path,  its  period,  control,  and  the  asso¬ 
ciated  performance  index.  In  the  last  section  of  this  chapter,  the 
analytical  results  are  compared  with  linear  extrapolations  near  the 
static  equilibrium  solution  and  with  numerical  results  of  the  previous 
chapter  in  the  region  near  the  local  minimum  operating  point. 

7. 1  Formulation  of  the  Problem 

The  perturbation  technique  which  is  used  in  this  chapter 
to  determine  an  approximate  analytical  solution  for  the  optimal 
periodic  control  problem,  equations  (5.1)  through  (5.5),  is  most 
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frequently  credited  to  Lindstedt  and  Poincare  [40,41].  Nineteenth 
century  astronomers,  such  as  Lindstedt  (1882),  Bohlin  (1889),  and 
Gylden  (1893),  developed  techniques  to  avoid  the  appearance  of  secular 
terms  in  perturbation  solutions  of  a  class  of  differential  equations. 
Poincare  (1892)  proved  that  the  expansions  obtained  by  Lindstedt 's 
technique  are  asymptotic.  In  this  section,  the  example  problem  is 
reformulated  in  a  similar  manner  by  expanding  the  independent  variable 
determining  a  functional  relationship  for  the  expansion  parameter,  and 
assuming  an  asymptotic  solution. 

To  simplify  the  perturbation  analysis,  the  two  point  bound¬ 
ary  value  problem,  equations  (5.7)  through  (5.10)  is  rewritten  in  the 
equivalent  form  of  a  single  fourth  order  differential  equation, 

.4  a 2  \2,2 

b — 4  =  *  — o  +  3(dt)  — 2  ”  y’  (7.1) 

dt  dr  V  '  dr 


where  the  new  variable  is  defined  by  y  s  x^.  The  periodicity  con¬ 
dition  is  applied  to  y  and  its  first  three  derivatives  at  the  bound¬ 
ary  points  corresponding  to  an  unrestricted  period,  T.  The  remaining 
variables  in  the  original  problem  may  be  obtained  from  the  solution 
to  equation  (7.1),  its  derivatives,  and  the  following  relationships 
derived  from  the  Euler-Lagrange  equations, 


(7.2) 


(7.3) 
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(7.4) 


The  general,  analytical  solution  to  equation  (7.1)  developed 
by  the  perturbation  technique  is  assumed  to  exist  in  the  form, 


y(e,t,ui)  =  y0(T,u>)  + 


e2y, 


(t,iu)  + 


(7.5) 


where  T  is  the  expanded  independent  variable  typically  defined  by  the 
expression 

t  =  t(1  +  eujj  +  e^u>2  +...),  (7.6) 

ui  is  the  infinite  set,  {ui^,  ui^,  . ..},  of  undetermined  constant  parame¬ 
ters  which  further  scale  the  independent  variable  in  (7.6),  and  £  is 
the  small,  perturbation  (expansion)  parameter.  The  general  solution, 
(7.5),  is  obtained  by  sequentially  determining  the  functional  represen¬ 
tation  of  its  components  related  to  successive  orders  of  the  expansion 
parameter.  In  this  manner,  the  approximation  may  be  developed  to  what¬ 
ever  order  is  desired. 

Before  proceeding  with  the  solution  development,  a  particular 
solution  to  equation  (7.1)  about  which  the  general  solution  is  to  be 
asymptotically  expanded  must  be  selected,  and  the  functional  relation¬ 
ship  of  the  expansion  parameter  to  the  system  parameters  must  also  be 
determined.  The  steady-state  solution  of  the  example  problem,  from 
the  results  in  section  5.2  of  the  frequency  test,  was  shown  to  be 
optimal  for  values  of  the  control  weighting  parameter,  b  >  <(.  For 
values  less  than  this  critical  value,  periodic  solutions  can  be  found 
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which  provide  better  performance.  With  the  objective  of  developing  a 
general  solution  for  the  example  problem  that  represents  the  class  of 
optimal  periodic  solutions,  an  appropriate  selection  for  the  particu¬ 
lar  solution  is  the  steady-state  solution  coupled  with  a  functional 
relationship  for  the  expansion  parameter  that  expresses  the  range  for 
b  <  fc.  The  following  expression  for  e  provides  that  relationship, 


t  -  1  "  e 

b  =  -5- 


(7.7) 


where  e  £  1.  The  quadratic  form  of  e  is  required  here,  as  will  become 
obvious  in  the  next  section,  to  permit  forming  expressions  that  annihi¬ 
late  secular  producing  terms  in  the  forcing  functions  of  the  expansion. 

Now  an  expression  for  the  differential  equation  (7.1)  may  be 
established  in  terms  of  the  expansion  parameter  e,  the  expanded  time 
variable  t,  and  the  infinite  set  of  constant,  undetermined  parameters 
u).  Since  T  is  linearly  related  to  t,  and  the  parameters  e  and  u>  are 
constants,  the  derivatives  of  y(e,T,u>)  with  respect  to  t  in  equation 
(7.1)  are 


&  _  dx 

dt  ~  dx  dt  ’ 


(7.8) 


d/d£  dx  \  d2y/dx \2 
'  dx\dx  dt  /  *  .  2\dt  /  » 


(7.9) 


dT4U/ 


(7.10) 


An  approximate,  analytical  solution  to  the  differential 


equation  (7.14),  which  is  equivalent  to  the  Euler-Lagrange  equations 


•  >  %  •  » 
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of  the  example  problem,  is  developed  in  this  section.  This  provides  a 
general  expression  for  periodic  extrema  of  the  optimal  periodic 
control  problem.  Then  the  components  of  the  solution  are  determined 
sequentially  by  substituting  the  assumed  solution,  (7.5),  into  the 
above  differential  equation,  (7.14),  and  equating  like  powers  of  the 
expansion  parameter,  e,  in  the  resulting  infinite  series  expansion. 
Operator  notation  for  the  derivative  with  respect  to  T,  i.e.,  D  =  d/dt, 
is  used  in  the  remainder  of  this  chapter  as  a  convenience  in  represent¬ 
ing  lengthy,  complex  differential  equations. 

By  equating  the  coefficients  of  in  the  expansion  of  equa¬ 
tion  (7.14),  the  following  differential  equation  for  the  yQ  component 
of  the  general  solution  is  obtained, 

D4yQ  =  -4D2y0  +  12(Dy0fy\  -  4y0.  (7.15) 

As  was  identified  in  the  previous  section,  a  solution  to  (7.15)  is  the 
steady-state  or  static  solution, 

y0  =  °yo  =  °2yo  =  °3yo  =  °-  (7-16) 

This  also  represents  the  particular  solution  about  which  the  general 
solution  is  expanded. 

Differential  equations  for  the  remaining  components  of  the 
general  solution  are  determined  in  the  same  manner  as  those  for  y^. 

Each  of  these  differential  equations  has  the  following  form, 

(D4  ♦  4D2+  4)y.  =  Fi(yi.1,...,y0), 


(7.17) 
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where  represents  the  forcing  function  of  the  i^  equation  in  terms 
of  previously  determined  components.  The  left  hand  side  of  the  equa¬ 
tion,  identical  for  each  unknown  component  y^,  identifies  the  charac¬ 
teristic  equation.  The  subscript,  i,  associated  with  the  unknown 
component,  y^,  and  with  the  forcing  function,  F^,  also  indicates  the 
power  of  the  expansion  parameter  whose  coefficients  in  the  infinite 
series  expansion  form  the  corresponding  differential  equation  (7.17). 
Forcing  functions  associated  with  the  unknown  components,  y  ,  through 
y^,  are  listed  below  after  simplifying  the  expressions  by  incorporating 
the  steady-state  solution  (7.16), 


F1  =  0,  (7.18) 

F2  =  4u>1D2(D2+2)y1,  (7.19) 

F3  =  4u)JD2(D2+2)y2  +  (l-10ui2+4u>2)D4yi 

-  (12u)2-8u>2)D2yi  +  12(Dy1)2D2y1,  (7.20) 

F4  =  4ui1D2(D2+2)y3  +  (l-10u>2+4u>2)D4y2 

-  (12u»2-8u>2)D2y2  +  12(Dy2)2D2y2 

♦  24Dy2Dy1D2y1  -  481^(1^  iVyj 

3  4 

+  (20u»1-20u)1u)2+4u>3-4u»1)D  yj 

+  (16uj^-24u)  w_+8u»9)D2y1 . 


(7.21) 
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Expressions  corresponding  to  any  desired  number  of  terms  in  the  gen¬ 
eral  solution  (7.5)  may  be  obtained  in  the  same  manner.  However,  the 
forcing  function  becomes  increasingly  more  complex  as  additional  terms 
are  added. 

To  obtain  expressions  for  the  components  of  the  general  solu¬ 
tion,  consider  first  the  solution  to  the  characteristic  equation, 

Vc  _  _ 

y.  =  (A.+  C.t)sinV2t  +  (B.+  Dix)cosV2t  (7.22) 

where  *  distinguishes  this  homogeneous  solution  from  the  complete 
solution.  In  order  to  satisfy  the  periodicity  conditions,  the 
secular  terms  in  (7.22)  must  vanish.  This  requires  that  the  arbitrary 
constants  and  must  be  chosen  identically  zero. 

Now  consider  the  contribution  to  the  component  solutions 
associated  with  the  forcing  terms.  Combining  this  particular  solution 
with  the  previous  homogeneous  result  (7.22)  determines  the  functional 
expression  for  respective  components  of  the  general  solution.  Since 
there  is  no  forcing  term  associated  with  the  equation  for  y^,  the 
expression  for  this  component  is 

yj  =  AjSin^t  +  BjCos^T.  (7.23) 

2 

The  forcing  function  (7.19)  is  also  zero  since  (D  *2)y^  =  0.  This 
yields  a  similar  result  for  the  y2  component 

y2  =  A2sinV2t  +  B2cosV2t  .  (7.24) 

Using  the  results  (7.23)  and  (7.24)  in  the  expression  (7.20)  for 
provides  the  following  relationship, 
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F3  =  -MAWj-l+aCA^+B^KAjSinVat  +  Bjcos^l) 

-laCAj-SB^AjSinSVat  -  UOA^-B^BjCos^T.  (7.25) 

Observe  that  the  first  term  in  the  forcing  function  is  a  solution  to 
the  characteristic  equation.  This  will  produce  a  secular  term  in  the 
solution  of  equation  (7.17).  Similar  secular  producing  terms  occur 
in  each  equation  for  successive  components  of  y.  To  satisfy  the  peri¬ 
odicity  condition,  the  arbitrary  constants  and  undetermined  parameters 
are  restricted  in  a  manner  that  resultant  coefficients  of  sia-j2x  and 
cosV2t  in  each  of  the  corresponding  forcing  functions  are  reduced  to 
zero.  The  restricting  equation  generated  by  the  forcing  function 
(7.25)  is 

-  1  ♦  3(Aj  +  Bj)  =  0.  (7.26) 

The  solution  for  y^  then  becomes 

=  AjiinV^T  ♦  B^cos^Zt  -  ^|(Aj-3Bj)A1sin3V2t 

♦^|(3a5-Bj)BiCos3V2t.  (7.27) 

Using  the  preceding  results  for  the  components  y^,  y2 ,  and  y^  in  the 
expression  (7.21)  for  determines  a  second  restricting  equation, 

3(AjA2  ♦  B^)  ♦  2ui1(w*  ♦  2 ui2)  =  0,  (7.28) 

and  a  solution  y^  to  equation  (7.17), 


y4  =  A^sinV2t  +  B^cos^X 

-  yfg  [(w1A1+6A2)(A^-B^)-2A1B1  (^8^682)18111341 

-  j|g( (WjBj^Bj) (A^-Bj)+2A1B1(w1A1+6A2)]cos34x.  (7*29) 

As  this  process  is  continued,  one  new  restricting  equation 
is  determined  for  each  additional  term  evaluated  in  the  expansion. 

It  should  be  noted  that  the  restricting  equation  which  relates  param¬ 
eters  through  the  i^*1  term  of  the  expansion  is  determined  from  the 
forcing  function  of  the  differential  equation  associated  with  the 
coefficients  in  the  expansion  of  the  i+2n<*  term. 

4 

An  approximate  solution  through  order  £  may  now  be  written 

in  the  form  of  equation  (7.5)  using  the  component  expressions  (7.23), 

(7.24),  (7.27),  and  (7.29).  The  resulting  expression  may  be  simplified 

with  essentially  no  loss  in  generality  by  imposing  the  arbitrary  bound¬ 
ary  condition  y(0)  =  0.  This  fixes  the  time  and  phase  relationship  of 

the  solution,  but  it  does  not  effect  its  amplitude  or  period.  As  a 

result  the  approximate  solution  becomes 

y(€,T,w)  =  eBjCos-^T  +  e  a ^cos^x 

+  e3[B3cosV2T  +  g|B3cos3V2t] 

+  £4[B4cos4t  +  ~(u>1B3+6BjB2)cos3V2T,  (7.30) 

4 

where  the  expanded  time  is  given  through  order  £  by 


t [  1  -  etUj  +  £2(u>^-u>2)  -  £3(uj^-2u>1u)2+uj3) 


+  e*(u»j-3u>2u>2+2w1u>3-ui£)  ] .  (7.31) 


By  rearranging  terms  in  equations  (7.30)  an  expression  in  the  form  of 


a  Fourier  series  is  obtained, 


y(e,T,u))  =  [eB^  +  e2B2  +  £3B3  +  e^B^]cosi/2t 


+  {£3B1  +  £S(w1B1+6B2)]^BjCos3V2t.  (7.32) 


Continuing  the  expansion  further  produces  additional  cosine 
terms  with  arguments  that  are  odd  integer  multiples  of  -J2X.  An  inter¬ 
esting  characteristic  of  the  expansion  is  that  the  smallest  power  of 
£  in  the  coefficient  of  cos  n-JlX  is  the  odd  integer  n.  This  property 
is  commonly  referred  to  as  the  D'Alembert  characteristic  [42],  and  it 
assures  the  expansion  satisfies  the  necessary  condition  for  conver¬ 
gence  of  the  series,  i.e.,  the  nt^1  term  approaches  zero  as  n  becomes 
increasingly  large. 

It  should  be  noted  at  this  point  that  the  approximate 
solution  to  order  en  has  2n  unknown  parameters,  B^  and  uj^,  where 
k  -  1,  2,  ...,  n.  However,  only  n  restricting  equations  relating 
these  parameters  can  be  determined  leaving  n  of  them  arbitrary.  In 
the  next  section,  n  additional  relationships  are  determined  by  apply¬ 
ing  the  conditions  for  an  optimal  period,  allowing  the  parameters  to 
be  evaluated. 


,  -  ^  i 


7 . 3  Optimal  Solution  with  Respect  to  Period 
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A  general  expression  was  derived  in  the  previous  section  for 
extrema  of  the  example  optimal  periodic  control  problem.  The  special 
transversality  condition  (5.16)  associated  with  an  unrestricted  period 
is  applied  in  this  section.  First  an  expansion  in  powers  of  £  is  de¬ 
rived  for  both  the  variational  Hamiltonian  and  the  performance  index. 
Then  the  expressions  are  set  equal  by  equating  coefficients  of  like 
powers  of  £.  This  results  in  n  new  relationships  for  the  parameters  of 
the  system.  Combined  with  the  n  restricting  equations  previously 
obtained,  all  2n  parameters  in  the  expression  for  extrema  may  be 
explicitly  determined.  As  a  result,  the  solution  (7.32)  represents  the 
first  n  terms  of  an  asymptotic  expansion  that  approximates  the  optimal 
solution  with  respect  to  an  unrestricted  period. 

The  variational  Hamiltonaian  and  the  performance  index  must 
be  evaluated  along  the  extremal  path  to  properly  apply  the  condition 
for  optimal  period.  The  extremal  path  expressed  in  terms  of  the 
original  problem  may  be  determined  from  the  definition  y  =  x^,  the 
extremal  solution  (7.33),  and  repeated  use  of  the  Euler-Lagrange 
equations,  (5.7)  through  (5.10).  This  results  in  the  following  series 
expansions  which  approximate  the  extremal  path, 


x 1  =  [eB1+£2B2+£3B3)cosV2T  +  e3g|B3cos3V2T,  (7.33) 

x2  =  -V2[£B1+£2(B2-w1B1)+£3(B3-io1B2+(u»2-u)2)B1)]sinV2T 

-£3^|V2B3sin3V2t, 


‘•W-  - 
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A.j  =  -|V2[eB1+e2(B2+w1B1)+e3(B3+tu1B2+u»2B1)]sinV2T 

■e3T28^2Blsin3^2t’  (7.35) 

X2  =  |leB1+e2(B2-2w1B1)+£3(B3-2w1B2+(3w2-2ui2-l)B1)]cosV2T 

+e3Yf|  B3cos3V2t,  (7.36) 

3 

where  terms  of  higher  order  than  e  have  been  truncated.  Using  the 
expression  for  the  optimal  control  (5.15),  the  functional  representa¬ 
tion  for  the  expansion  parameter  (7.7),  and  the  expression  (7.36) 
provides  the  series  representation  of  the  extrema  control 

u  =  -2  [eBl+£2(B2-2U)1B1)+e3  (82-2^2+ (3u>2-2io2-l)B1)]cosV2x 

-e3  ||  b|cos3V2i.  (7.37) 

Substituting  these  expressions  in  equation  (5.6),  the  follow¬ 

ing  general  expression  for  the  variational  Hamiltonian  evaluated  along 
the  extremal  path  is  obtained, 

H  =  c32w1b|  +  e4[4ui1B1B2+(2u>2-3ui|)b|+  |b4] 

*£5 {2u)j (b|+2B^B3~  |b|)-(9ui|-4u>2-  |)BjB2 

+(7u»|-6u)1U)2+ui3)b|]  .  (7.38) 


•  »  .»•  •  •• • 


The  first  three  restricting  equations  were  used  to  obtain  this  simpli- 

4 

fied  form  of  the  Hamiltonian.  Only  terms  upto  £  in  the  expressions 

(7.33)  through  (1.31)  were  required  to  obtain  the  term  in  (7.38). 
Note,  too,  that  the  functional  representation  of  the  variational 
Hamiltonian  is  a  constant  expression,  as  it  should  be,  when  evaluated 
along  an  extremal  path. 

In  a  similar  manner,  a  general  expression  for  the  performance 
index  evaluated  along  the  extremal  path  is  determined  by  substituting 

(7.33) ,  (7.34),  and  (7.37)  in  the  expression  for  the  performance  index 
(5.1)  and  using  the  orthogonality  characteristics  [43]  of  the  sine  and 
cosine  in  performing  the  indicated  integration.  This  results  in  the 
following  constant  expression, 


J  =  lB?  -£5  fB?(B2-u>1B1). 


(7.39) 


It  is  interesting  that  the  improvement  in  the  performance  is  zero  to 
the  fourth  order. 

With  expressions  for  the  Hamiltonian  and  the  performance 
index  both  evaluated  along  an  extremal  path,  the  condition  for  optimal 
period  may  be  applied.  Equating  the  coefficients  of  like  powers  of  £ 
in  H  and  J  provides  the  additional  equations  needed  to  specify  the 
remaining  arbitrary  parameters  in  the  solution.  The  first  three 
of  these  optimal  period  relationships  are: 


0, 


(7.40) 
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tojBjBj  +  (2u)2-3u)^)B^+  |bA  =  0,  (7.41) 

2u)j (B2+2B1B3“B j)  -  (11^-4u)2-2)B1B2 

+  (910^-6^2 +u»3)B^  =  0.  (7.42) 


The  first  three  restricting  equations,  (7.26),  (7.29),  and  that 
from  the  expansion  to  are  rewritten  below,  using  the  results, 

Ai  =  0  for  i  =  1,  2,  ...,  obtained  in  the  previous  section  by  impos¬ 
ing  the  boundary  condition  y(0)  =  0.  The  simplified  restricting 
equations  are: 


4uiJ  -  1  +  3Bj  =  0, 

38^2  +  2u>1(uy»-2u)2)  =  0, 

97  S  2  2 

Wl  '  2AB?3  "  12B1B2 

2  2  2 
-  4B1(w1+12u)1u>2+8w1u>3+4u>2)  =  0. 


(7.43) 

(7.44) 


(7.45) 


The  six  equations,  (7.40)  through  (7.45),  may  be  solved  to 
obtain  values  for  the  constant  parameters  that  determine  explicitly 

3 

through  e  the  approximate  analytical  solution  that  satifies  all 
the  first  order  necessary  conditions  for  an  optimum.  Equation  (7.40) 
allows  the  choice  of  either  B^  =  0  or  tUj  =  0.  The  initial  choice 
leads  to  Bj  =  B2  =  B3  =  0,  which  results  in  the  static  equilibrium 
solution.  For  z  <  k,,  the  static  solution  is  a  local  maximum.  For 
the  second  choice,  u)j  =  0,  the  parameters  are  determined  to  be 
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1 

4* 

(7.46) 

Note  that  the  parameter  has  two  values.  As  more  terms  are 
included  in  the  expansion,  the  additional  restricting  equations 
and  the  optimal  period  relationships  become  increasingly  complex. 

This  results  in  multiple  valued  parameters  which  indicate  the 
possibility  of  obtaining  analytical  expressions  for  numerous 
locally  optimal  solutions. 

Using  the  values  for  the  parameters  listed  in  (7.46)  leads 
to  the  following  expressions  in  terms  of  e  and  t  for  the  optimal 
path  and  control. 


j  =  0,  Bj  =  ±  ^|,  B2  =  0,  U>2 


u.3  =  0,  and  B3  -  .  m 


= ;  hd 

1  1  AA  • 


X1  =  ±[  (^3e  ■  |^|£3)cosV2t  +  ^^e^cosa^T] ,  (7.47) 
X2  "  +I(^fe  +  ^|g|e3)sillV2T  +  ^|e3sin3V2T] ,  (7.48) 
\  =  ~  ^|g|e3)sinV2T  +  y^|2e3sin3V2T] ,  (7-49) 

K2  =  ±I(4e  '  ^g|£3)cosV2T  +  |^e3cos3V2l],  (7.50) 
u  s  +((^3^2  +  7273e3^COSl^T  +  J »  (7-51) 
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2 

where  the  expanded  time  may  be  expressed  ast=(l+fce)t.  The  above 
solutions  are  an  approximate  analytical  representation  o£  the  local 
minima  associated  with  the  principal  family  of  periodic  solutions 
identified  in  the  previous  chapter.  In  each  of  these  expressions 
terms  of  £  and  higher  have  been  truncated. 

In  a  similar  manner,  the  constant  expressions  in  terms  of  £, 
for  the  period,  variational  Hamiltonian,  and  performance  index  that 
correspond  to  the  previous  minimum  solution  are  found  to  be 


nj2(l 

12..  . 

-  ^r£  )  ,  and 

(7.52) 

J  = 

-  -V 

24e 

(7.53) 

Terms  of  order  £^  and  higher  have  been  truncated  in  the  derivation  of 
the  last  expression. 

7.4  Verification  of  Results 

The  anlytical  results  developed  in  this  chapter  are  verified 
in  this  section.  The  functional  representation  of  extremal  solutions 
is  compared  qualitatively  with  the  numerical  results  obtained  in  the 
previous  chapter.  In  the  region  near  the  static  equilibrium  point  the 
extremal  solution  is  compared  to  the  result  of  a  linear  analysis  about 
the  same  point.  Finally,  the  analytical  results  obtained  in  the  last 
section  for  the  optimal  solution  are  compared  with  corresponding  numer¬ 
ical  results  of  the  previous  chapter. 
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It  is  interesting  that  except  for  the  static  equilibrium 
solution  and  the  the  optimal  solutions  with  respect  to  an  unrestricted 
period,  the  constants  of  the  periodic  extremal  solutions,  given  by 
equations  (7.33)  through  (7.37),  are  indeterminant.  Examining  the 
extremal  solution  derived  from  an  expansion  up  to  eQ  shows  that  there 
are  2n  constants  to  be  determined.  However,  there  are  only  n  restrict¬ 
ing  equations  that  provide  relationships  between  these  constants.  As 
a  result,  n  of  them  may  be  arbitrarily  chosen.  This  provides  for 
large  regions  of  the  phase  space  to  be  densely  packed  with  periodic 
extrema  as  was  shown  to  be  the  case  for  the  computational  results 
obtained  in  the  previous  chapter.  Considering  additional  higher  order 
terms  in  the  asymptotic  expansion  greatly  complicates  the  associated 
relationships  for  determining  the  parameters.  As  the  restricting 
equations  and  those  for  optimal  period  become  more  complex,  the  new 
constants,  determined  by  these  equations  become  multi-valued.  This 
in  turn  leads  to  the  identification  of  multiple  local  extrema,  each  of 
which  satisfies  the  optimal  period  condition.  Results  obtained  numer¬ 
ically  in  the  previous  chapter  are  in  agreement. 

Now  consider  the  linear  region  near  the  static  equilibrium 
point.  Choosing  the  parameters  !$£  =  =  0  and  allowing  ■*  0  in 

the  expression  for  the  solution  restricts  it  to  the  region  near  the 
static  path.  Using  these  values  in  the  restricting  equations,  (7.43) 
through  (7.45),  determines  the  remaining  parameters, 

W1  =  1  b  w2  =  '  b  w3  *  ±  l£*  (7'54) 
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The  fundamental  frequencies  defined  by  the  solution  given  by  equa- 

3 

tions,  (7.31)  and  (7.32)  can  be  written  through  the  £  term  as 

u>  =  V2(l  ±  +  |e2  ±  y|e3).  (7.55) 


Precisely  the  same  result  is  obtained  for  the  fast  and  slow 
frequencies  predicted  by  the  eigenvalues  of  the  linearized  Euler- 
Lagrange  equations.  The  coefficient  matrix  for  this  set  of  linear 
equations  is 


■  0  1  0  0" 
0  0  0  -i 

-10  0  0 
.0  1-1  0. 


(7.56) 


evaluated  at  the  static  equilibrium  point.  The  eigenvalues  of  the 
coefficient  matrix  are  given  by 


M2  =  2S  l±>/1-4b  '  (7.57) 

Using  the  relationship  between  the  control  weighting  parameter,  b,  and 
the  expansion  coefficient,  e,  then  simplifying  the  result  gives  the 
following  expression  for  the  frequency, 

u>  =  V2(l  ±  e)‘V  (7.58) 

where  p  =  iu>.  Expanding  this  expression  using  the  generalized  bi¬ 
nomial  theorem  gives  precise  agreement  with  equation  (7.55). 


Finally,  a  comparison  is  made  in  table  7.1  between  the  asymp¬ 


totic  approximation  of  the  optimal  solution  and  the  corresponding 
computational  result  for  several  values  of  the  expansion  parameter. 

In  the  table,  x} ,  and  X^  identify  the  initial  values  of  the  optimal 
periodic  solution,  where  both  and  X^  have  initial  value  zero.  The 
optimal  period  is  designated  T  and  the  performance  index,  J. 

The  results  of  the  asymptotic  expansion  after  three  terms 
provides  excellent  agreement  with  the  computational  results  for  the 
smaller  values  of  the  expansion  paramter  e.  It  provides  a  good  first 
approximation  for  large  values  of  e.  The  agreement  obtained  is  much 
better  than  would  be  predicted  by  assuming  an  accuracy  equivalent  to 
the  order  of  e  in  the  first  truncated  term. 
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Table  7.1  Comparative  Results 


CHAPTER  8 


CONCLUSIONS  AND  RECOMMENDATIONS 

An  attempt  has  been  made  in  this  study  of  optimal  periodic 
control  theory  to  merge  pertinent  experience  and  theory  from  the  fields 
of  analytical  dynamics  and  celestial  mechanics  with  that  existing  in 
the  optimization  and  control  communities.  This  has  resulted  in  some 
important  theoretical  contributions  for  this  special  class  of  optimal 
control  problems.  It  has  also  developed  a  general  characterization  of 
optimal  periodic  control  processes  evolving  from  static  equilibrium 
solutions  of  this  problem.  A  summary  of  the  more  significant  results 
and  recommendations  for  future  study  are  presented  in  this  chapter. 

8. 1  Summary  of  Conclusions 

The  principal  contributions  to  the  general  theory  of  optimal 
periodic  control  are  summarized,  and  particular  results  of  both  a  nu¬ 
merical  and  analytical  investigation  of  an  optimal  periodic  control 
problem  are  reviewed  in  this  section.  There  are  two  key  relationships 
developed  in  this  study  which  are  basic  to  most  of  the  theoretical 
results.  One  is  an  algebraic  equation  which  relates  the  initial 
conditions  for  periodic  solutions  to  the  Riccati  differential  equation 
with  the  elements  of  the  monodromy  matrix.  The  other  is  a  particular 
similarity  transformation  of  the  monodromy  siatrix  which  results  in  a 
partitioning  of  its  eigenvalues  into  the  two  submatrices  on  its 
principal  diagonal. 
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A  relationship  between  real  values  of  the  Riccati  variable 
and  the  eignevalues  of  the  aoaodromy  matrix  is  derived  by  exploiting 
the  symplectic  properties  of  the  monodromy  matrix,  and  the  special 
properties  of  its  eigenvalues  both  preserved  through  the  similarity 
transformation.  The  important  conclusion  derived  from  this  relation¬ 
ship  is  that  locally  optimal  periodic  solutions  of  the  control  problem 
are  generally  unstable.  If  an  eigenvalue  associated  with  an  optimal 
periodic  solution  lies  on  the  unit  circle  at  other  than  1  or  -1,  it 
must  occur  in  sets  of  double  conjugate  pairs.  At  the  critical  points 
1  and  -1,  the  eigenvalues  must  occur  in  coupled  pairs. 

A  deficiency  is  corrected  in  the  previously  existing  suffi¬ 
cient  conditions  for  local  optimality  of  periodic  control  problems 
with  fixed  periods.  The  specific  correction  identifies  that  two  unity 
eigenvalues  are  associated  with  all  extrema  of  the  optimal  periodic 
control  problem.  This  is  a  property  of  the  monodromy  matrix  of  a 
Hamiltonian  system  such  as  the  Euler-Lagrange  equations  in  the  optimal 
periodic  control  problem.  The  sufficiency  conditions  are  then  extend¬ 
ed  to  include  the  class  of  problems  for  which  the  period  is  unrestrict¬ 
ed.  The  neighboring  optimal  control  law,  derived  in  the  above  process, 
is  used  to  develop  a  periodic  regulator.  It  optimizes  the  return  path 
for  small  perturbations  from  the  nominal  path  in  a  manner  analagous  to 
the  static  regulator. 

Two  very  interesting  and  simple  results  are  derived  from  the 
preceding  general  developments  for  symmetric,  fourth  order  optimal 
control  problems.  An  inequality  expression  involving  only  the  trace 
of  the  monodromy  matrix  provides  a  necessary  condition  for  optimal 
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periodic  solutions.  The  second  result  provides  an  explicit  expression 
in  terms  of  elements  of  the  monodromy  matrix  for  the  initial  conditions 
of  periodic  solutions  to  the  Riccati  differential  equation.  This  is 
used  to  show  that  only  two  periodic  solutions  of  the  Riccati  equation 
can  exist  for  such  a  system. 

In  the  second  half  of  this  study  a  particular  optimal  control 
problem,  constructed  to  characterize  optimal  periodic  solutions,  is  in¬ 
vestigated  both  numerically  and  analytically.  The  remarkable  complex¬ 
ity  of  its  solutions  is  characterized  by  the  numerical  results.  An 
infinity  of  extremal  solutions  are  found  forming  continuous,  one- 
parameter  families  that  densely  pack  large  regions  of  the  phase  space. 
They  intersect  at  common  solutions  which  may  be  determined  by  a  measure 
of  their  stability.  Characteristics,  such  as  the  number  of  arcs  in  a 
period,  are  identified  to  distinguish  solutions  of  different  families. 
The  computer  program  used  in  this  study  exploits  these  distinguishing 
characteristics  and  obtains  extremal  solutions  by  tracing  families  in 
a  systematic  and  predictable  manner.  Multiple  solutions  are  found  that 
satisfy  the  necessary  condition  for  an  optimal  period.  The  Riccati 
variable  associated  with  these  solutions  is  found  to  exist  over  the 
full  period  for  only  the  solution  associated  with  the  principal  family. 
However,  it  is  conjectured  that  other  solutions  exist  that  also  satisfy 
each  condition  for  local  optimality.  Similar  results  for  local  optima 
with  specified  period  are  determined. 

The  optimal  periodic  control  problem  is  also  analyzed  using 
a  perturbation  technique  which  results  in  reformulating  the  problem  as 
an  asymptotic  expansion.  General  expressions  for  the  solution's  state, 


period  and  performance  index  is  derived  in  terms  of  the  independent 
variable,  time,  and  an  expansion  parameter.  The  expressions  are  in 
the  form  of  a  Fourier  series  with  its  fundamental  frequency  and  coef¬ 
ficients  represented  by  asymptotic  expansions.  A  comparison  of  the 
analytical  results  with  relatively  few  terms  to  the  numerical  results 
is  quite  good. 

8.2  Recommendations  for  Future  Study 

There  are  several  avenues  along  which  future  effort  might  be 
directed.  The  principal  problem  that  motivated  this  study  was  the 
apparent  lack  of  convergence  of  typical  first  order  numerical  optimiza¬ 
tion  techniques  for  particular  optimal  periodic  control  problems.  The 
development  of  efficient  computational  algorithms  for  this  class  of 
optimization  problem  is  an  important  and  fruitful  path.  Higher  order, 
more  complex  systems  must  be  addressed  such  as  in  the  recent  work  on  a 
hypersonic  cruise  cycling  trajectory,  by  D.  Walker  [44]. 

The  numerical  investigation  initiated  in  this  study  is 
virtually  limitless.  Solutions  with  initial  conditions  lying  in  only 
one  region  of  the  initial  condition  surface  of  symmetry  was  investi¬ 
gated  in  this  study.  Nonsymmetric  solutions  exist  and  theory  is  avail¬ 
able  to  predict  intersections  with  symmetric  solutions.  Also,  irregu¬ 
lar  families  that  have  no  interconnection  at  all  with  the  principal 
families  have  been  shown  to  exist  for  other  dynamic  problems  [25]. 
Additional  optimal  solutions  with  greater  improvement  in  performance 
may  be  found.  Another  possible  avenue  of  research  is  the  investigation 
of  quasi-periodic  solutions  which  may  provide  better  performance  for 
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the  infinite  time  problem  (not  addressed  in  this  study)  which  is  a 
generalization  of  a  periodic  solution  of  increasingly  large  period. 

The  further  development  of  the  asymptotic  series  expansions 
should  identify  additional  optimal  solutions.  It  might  also  lead  to 
the  generalization  of  the  frequency  test  to  determine  optimality  for 
arbitrary  extremal  trajectories. 
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APPENDIX  A 


A  SOLUTION  TO  THE  RICCATI  DIFFERENTIAL  EQUATION 


By  solving  an  appropriate  set  of  linear  differential 
equations,  a  solution  to  the  following  Riccati  differential  equation 
may  be  obtained, 

P(t)  =  -AT(t)P(t)-P(t)A(t)-P(t)B(t)P(t)+C(t).  (A. I) 


The  following  is  essentially  a  restatement  of  a  Lemma  and  proof 
presented  by  Brockett  [36],  pages  155-157. 

Given  the  set  of  2n  linear  differential  equations, 


aJ 


A(t)  B(t) 


C(t)  -A  (t)JLAj 


(A. 2) 


let  its  transition  matrix,  ♦,  be  partitioned  into  four  n-square  sub¬ 
matrices  as  follows, 


*(t,t0)  = 


(b  <b 
11  v12 


<t>  <t> 

.21  22. 


(A. 3) 


Then  a  solution  to  equation  (A.l)  is 

P(t)  =  [*21  +  *22P(t0)][4>n  +  <*,12P(t0)l~1  (A. 4) 

when  the  indicated  inverse  exists. 
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In  order  to  show  (A. 4)  is  in  fact  a  solution  to  (A.l),  first  examine 
the  boundary  conditions.  Applying  4>(tQ,tQ)  =  I  to  equation  (A. 4) 
gives  P(tp)  =  P(tg).  Clearly  the  boundary  conditions  are  satisfied. 
Substituting  (A. 4)  into  (A.l)  and  reducing  the  result  to  an  identity 
would  complete  the  verification.  To  simplify  the  algebra  involved  in 
doing  this,  first  define  the  n-square  matrices  X  and  A  by 


x(t) 

A(t) 


(A. 5) 


Since  41  is  a  solution  to  equation  (A. 2),  the  linear  combination  of 
4>  in  (A. 5)  is  also;  hence, 


”xct) 

A(t) 

B(t)  " 

X(t) 

A(t)_ 

C(t) 

w 

-AT(t)_ 

A(t). 

Using  X  and  A  from  (A. 5),  the  expression  for  the  solution  P(t)  given 
by  (A. 4)  can  be  written, 

P(t)  =  A(t)X_1(t).  (A. 7) 


Since  the  time  derivative  of  the  inverse  of  a  matrix  may  be  expressed 


X_1(t)  =  -x'1(t)X(t)X-1(t), 


(A.  8) 


the  time  derivative  of  equation  (A. 7)  is 
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P(t)  =  A(t)X_1(t)  -  A(t)X_1(t)X(t)X'1(t)  .  (A. 9) 

Using  the  expressions  for  X(t)  and  A(t)  from  (A. 6),  and  for  P(t)  from 
(A. 4),  the  expression  in  (A. 9)  is  reduced  to  the  expression  in  (A.l) 
as  follows, 

P(t)  =  C(t)X(t)X*1(t)  -  AT(t)A(t)X_1(t) 

-A(t)X_1(t)A(t)X(t)X_1(t)  (A. 10) 

-A(t)X_1(t)B(t)A(t)x'1(t), 

P(t)  =  C(t)  -  AT(t)P(t)  -  P(t)A(t)  -P(t)BP(t) .  (A. 11) 

Thus  it  has  been  shown  that  the  expression  in  (A. 4)  is  the  general 
solution  to  the  Riccati  equation  (A.l). 


APPENDIX  B 


THE  MONODROMY  MATRIX  OF  A  SYMMETRIC  SYSTEM 

The  following  derivation  of  the  monodromy  matrix  for  a 
symmetric,  fourth  order  system  is  an  extraction  from  the  notes  of 
Dr.  Roger  A.  Broucke. 

Assume  a  fourth  order,  dynamical  system  of  the  form 


x  =  Ux(x,y),  x(t+T)  =  x(t), 

y  =  uy(x,y),  y(t+T)  =  y(t) ,  (B.l) 


where  x  and  y  are  scalars  and  U(x,y)  is  a  potential  function. 

Assume  further  that  the  system  is  symmetric,  satisfying 

x(-t)  =  x(t) ;  x(-t)  =  -x(t), 

y(-t)  =  -y(t) ;  y(-t)  =  y(t) ,  (B.2) 


such  that  the  system  is  invariant.  This  implies  the  following  rela¬ 
tionship, 

U(x,y)  =  U  (x,y)x(t)  +  U  (x,y)y(t)  =  0.  (B.3) 

a  y 

Using  equations  (B.l)  and  (B.2)  in  equation  (B.3)  results  in  the  fol¬ 
lowing  useful  expressions, 


ux(x,y)  =  ux(x,-y), 

U  (x,y)  =  -U  (x,-y) . 


(B.4) 
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The  variational  equations  of  (B.l)  can  be  written  as 


6x  =  U  fix  +  U  6y, 
xx  xy  ' 


6y  =  U  fix  +  U  6y, 
yx  yy 


which  applying  (B.l)  and  (B.2)  again  require  that 


(B.5) 


•-  1 


Uxx(x*y) 

=  Uxx(x*-y)* 

1  Uxy(x'y) 

=  -uxy(x,-y), 

Uyx(X>y) 

*  -Uyx(X,'y)> 

1  uyy(*’y) 

*  Uyy(X,'y)‘ 

(B.6) 

Rewriting  equations  (B.5)  in  the  form  of  four  first  order 
equations,  where 

vj  =  6x, 

v2  =  5y» 

v3  =  6x’ 

v4  *  6y- 

This  results  in  the  following  set  of  equations 

*1  =  v 

V,  =  V 


(B.7) 


V 

4  ’ 

1  V, 

+  U  v„ 

XX  1 

xy  2 

1  V, 

+  U  v„ 

yx  1 

yy  2 

(B.8) 


■L~y- '  ‘■•vv'1-'. 
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which  in  vector  notation  is 


V  =  A  V, 


(B.9) 


where  V  is  the  vector  composed  of  v^,  v^,  v^,  and  v^;  and  A  is  defined 


as 


A  = 


0  0  10 
0  0  0  1 


U  U  0  0 

xx  xy 

,  U  U  0  0 

L  y*  yy 


(B. 10) 


Using  the  symmetry  properties  (B.6),  it  is  easy  to  verify  that 

A(-t)  =  -BA(t)B,  (B. 11) 

where 


B  = 


10  0  0 

0-100 
0  0-10 
0  0  0  1 

Note  that  B  has  the  following  properties: 


(B. 12) 


T  T  -1 

BB  =  I;  B  =  B  =  B. 


(B. 13) 


Consider  the  transition  matrix  4(t,0)  for  this  system.  It 
satisfies  the  variational  equations  (B.8) 

♦(t,0)  =  A(t)4»(t,0).  (B.  14) 


Replacing  t  by  -t  and,  using  the  previous  symmetry  relations,  equation 
(B.14)  becomes 


4>(-t,0)  =  -A(-t)4>(-t,0), 


which  using  equation  (B.ll)  can  be  written 


(B. 15) 


4>(-t,0)  =  BA(t)B4>(-t,0)  or 


B4>(-t,0)  =  A(t)B<K-t,0). 


(B. 16) 


This  shows  B$(-t,0)  is  a  fundamental  solution  matrix  of  equation  (B.8) 
and  must  be  related  to  the  transition  matrix  (or  to  any  other  funda¬ 
mental  solution  matrix)  by  post  multiplication  of  a  constant  matrix, 


B4>(-t  ,0)  =  <J>(t,0)C. 


(B. 17) 


From  the  initial  conditions  the  constant  matrix  C  must  be  B,  therefore 


the  following  relationship  may  be  written, 


<K-t,0)  =  B«Kt,0)B. 


(B. 18) 


Finally,  these  symmetry  properties  applied  at  a  particular  value 
of  t,  equal  to  one  full  period  T,  are  considered.  From  properties  of 
the  mondromy  matrix,  the  following  relation  may  be  written, 


<Kt+T,0)  =  4>(t,0)*(T,0). 


When  t  =  -T,  equation  (B.19)  reduces  to 


(B. 19) 


I  =  4>(-T,0)<KT,0). 


From  equations  (B.20)  and  (B.18),  the  following  may  be  written, 


♦  (T,0)  =  *(-T,0)  =  B4>(T,0)B. 


(B.20) 


(B.21) 


"  '■'W' 
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Since  the  monodromy  matrix  is  symplectic, 

4>_1(T,0)  =  -K4>T(T,0)K,  (B.22) 

where  K  is  the  fundamental  symplectic  matrix.  Combining  equations 
(B.21)  and  (B.22)  gives 

-K4>T(T,0)K  =  B*(T,0)B, 

<t>T(T,0)  =  -KB4>(T,0)BK, 

4>T(T,0)  =  (BK)*(T,0)(BK),  (B.23) 


where,  from  the  definitions  of  B  and  K, 

”o  0  1  0~ 

0  0  0  >1 

BK  = 

10  0  0 

0-100 


(B.24) 


Forming  the  identities  which  result  from  equating  the  sixteen  elements 
of  the  resulting  two  matrices  from  equation  (B.23)  establishes  six 
nontrivial  symmetry  relationships  for  the  monodromy  matrix,  ♦(T,0), 


'33 

*11’ 

'*34  = 

*21’ 

'43  = 

*12’ 

11 

*22’ 

(B.25) 

'32 

*41’ 

11 

H 

0 

*23' 

From  these  relationships,  the  monodromy  matrix  for  a  symmetric,  fourth 
order  system  may  be  written  in  the  following  form, 
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♦Cr,o)  = 


a 

c 

h 

-i 


b  e  f 

d  -f  g 

i  a  -c 

j  -b  d 


(B.26) 


This  result  also  shows  that  only  the  first  two  columns 
of  the  monodromy  matrix  need  be  computed  to  obtain  the  trace  of  the 
matrix. 


APPENDIX  C 


CLASSIFICATION  OF  CRITICAL  SOLUTIONS 

Critical  points  (solutions)  of  the  first  and  second  kinds 
are  identified  and  classified  in  this  section  for  the  principal 
family  of  solutions  determined  in  the  numerical  investigation  of 
chapter  six.  This  is  directly  attributable  to  the  work  of  Henon 
[27]  and  Contopoulos  [26],  studying  the  stability  characteristics 
of  fourth  order  periodic  systems.  No  attempt  is  made  to  develop  or 
reproduce  the  associated  theory;  only  the  classification  by  type  of 
the  critical  solutions  of  the  sample  problem  and  the  significance  are 
expressed. 

On  the  principal  family ,  there  are  four  critical  points 
of  the  first  kind  (kc  =  +2);  1A,  1B1,  1B2,  and  1B3.  Three  critical 
points  of  the  second  kind  (kc  =  -2);  2A,  2B1,  and  2B2  also  exist. 

The  location  of  each  of  these  points  is  shown  in  figure  6.3  on  the 
initial  condition  plot  of  the  principal  family.  In  table  C.l  the 
initial  conditions  associated  with  the  critical  solutions  are 
identified.  The  Henon  matrix  and  resulting  classification  type  are 
specified  as  well.  The  Henon  matrix  is  essentially  the  result  of 
reducing  the  order  of  the  system  by  two.  Corresponding  to  the 
off-diagonal  elements  of  this  matrix,  different  types  of  exchange 
of  stability  are  identifiable.  For  the  critical  points  of  the  first 
kind,  each  is  associated  with  a  particular  Henon  type.  The  critical 

157 


59 


points  of  the  second  kind  are  all  of  the  type  not  classified  by  Henon 
(both  off-diagonal  elements  are  zero).  Contopoulos'  classification 
type  is  identified  by  these  cases. 

The  significance  of  the  classification  types  is  illustrated 
in  figure  C.l  for  critical  solutions  if  both  first  and  second  kinds. 

The  solid  line  identifies  unstable  solutions;  the  dashed  line  indicates 
the  stable  solutions;  and  the  dotted  line  represents  nonsynmetric 
solutions  that  do  not  lie  in  the  symmetric  initial  condition  surface. 
Each  critical  solution  is  found  at  the  intersection  of  families  of 
solutions  or  at  the  stability  transition  point. 

Solutions  1A  and  1B3  are  extremum  in  energy  (Hamilitonian) 
and  represent  only  the  transition  between  stable  and  unstable  solutions 
on  the  principal  family.  No  bifurcation  occurs  at  this  point  as 
depicted  in  figure  C.la.  Solution  1B1  is  a  bifurcation  point  and  is 
represented  in  figure  C.lb.  The  final  critical  solution  of  the  first 
kind,  1B2,  is  shown  in  figure  C.lc.  The  intersection  of  a  nonsynmetric 
family  with  the  principal  family  is  predicted  at  this  point. 

The  critical  solutions  of  the  second  kind  are  all  of  the 
type  D' ,  as  classified  by  Contopoulos.  The  upper  half  solutions 
(H  >  0)  were  not  investigated  and,  therefore,  are  not  further  classi¬ 
fied.  The  solution  2A,  shown  in  figure  C.ld,  is  identified  as  a  tri¬ 
furcation  point.  The  two  symmetric  families  have  been  traced.  The 
nonsynmetric  family  of  solutions  has  not  yet  been  found. 
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